{'fraction_non_alphanumeric': 0.045304172027200884, 'fraction_numerical': 0.017391104576364638, 'mean_word_length': 4.2130538922155685, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Combining Forecasts under Structural Breaks Using Graphical LASSO

In this paper we develop a novel method of combining many forecasts based on a machine learning algorithm called Graphical LASSO. We visualize forecast errors from different forecasters as a network of interacting entities and generalize network inference in the presence of common factor structure and structural breaks. First, we note that forecasters often use common information and hence make common mistakes, which makes the forecast errors exhibit common factor structures. We propose the Factor Graphical LASSO (Factor GLASSO), which separates common forecast errors from the idiosyncratic errors and exploits sparsity of the precision matrix of the latter. Second, since the network of experts changes over time as a response to unstable environments such as recessions, it is unreasonable to assume constant forecast combination weights. Hence, we propose Regime-Dependent Factor Graphical LASSO (RD-Factor GLASSO) and develop its scalable implementation using the Alternating Direction Method of Multipliers (ADMM) to estimate regime-dependent forecast combination weights. The empirical application to forecasting macroeconomic series using the data of the European Central Bank's Survey of Professional Forecasters (ECB SPF) demonstrates superior performance of a combined forecast using Factor GLASSO and RD-Factor GLASSO.


Introduction
A search for the best forecast combination has been an important on-going research question in economics. Clemen (1989) pointed out that combining forecasts is "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 Diebold and Shin (2019), 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, Diebold and Shin (2019) develop "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 (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 Bates and Granger (1969), 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 Graphical LASSO (Friedman et al. (2008)) and nodewise regression (Meinshausen and Bühlmann (2006)). 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. We show that graphical models fail to recover the entries of a nonsparse 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 latent factors. We call our algorithm the Factor Graphical LASSO. We use a factor model to estimate an idiosyncratic component of the forecast errors, and then apply a Graphical LASSO for the estimation of the precision matrix of the idiosyncratic component.
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 develops a unified framework to generalize network inference in the presence of structural breaks. 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 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, 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, 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 Regime-Dependent Factor Graphical LASSO (RD-Factor GLASSO) and develop its scalable implementation using ADMM to estimate regime-dependent forecast combination weights. Fourth, 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.
The paper is structured as follows. Section 2 reviews Graphical LASSO. Section 3 studies the approximate factor model for the forecast errors. Section 4 introduces Factor Graphical LASSO and contains theoretical results on the consistency of the Factor GLASSO estimator.
Section 5 introduces Regime-Dependent graphical model and discusses its implementation using ADMM. Section 6 validates theoretical results using simulations. Section 7 studies an empirical application for macroeconomic time-series forecasting. Section 8 concludes.
Notation. For the convenience of the reader, we summarize the notation to be used throughout the paper. Let S p denote the set of all p × p symmetric matrices. For any matrix C, its (i, j)-th element is denoted as c ij . Given a vector u ∈ R d and a parameter a ∈ [1, ∞), let u a denote a -norm. Given a matrix U ∈ S p , let λ max (U) ≡ λ 1 (U) ≥ λ 2 (U) ≥ . . . ≥ λ min (U) ≡ λ p (U) be the eigenvalues of U. Given a matrix U ∈ R p×p and parameters a, b ∈ [1, ∞), let |||U||| a,b ≡ max y a =1 Uy b denote the induced matrix-operator norm. The special cases are |||U||| 1 ≡ max 1≤j≤p p i=1 |u ij | for the 1 / 1 -operator norm; the operator norm ( 2 -matrix norm) |||U||| 2 2 ≡ λ max (UU ) is equal to the maximal singular value of U. Finally, U max ≡ max i,j |u ij | denotes the element-wise maximum.

Graphical Models for Forecast Errors
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, . . . , T . Let e t = (e 1t , . . . , e pt ) ∼ N (0, Σ) be a p × 1 vector of forecast errors. Assume they follow a Gaussian distribution. The precision matrix Σ −1 ≡ Θ contains information about partial covariances between the variables. For instance, if θ 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 W be the estimate of Σ. Given a sample {e t } T t=1 , let S = (1/T ) T t=1 (e t )(e t ) denote the sample covariance matrix, which can be used as a choice for W. Also, let Γ 2 ≡ diag(W) and its (i, j)-th element is denoted as γ ij . We can write down truncated Gaussian loglikelihood (up to constants) l(Θ) = log det(Θ) − trace(WΘ). When W = S, the maximum likelihood estimator of Θ is Θ = S −1 . The objective function associated with truncated Gaussian log-likelihood is also known as Bregman divergence and was shown to be applicable for non-Gaussian distributions (Ravikumar et al. (2011)).
In the high-dimensional settings it is necessary to regularize the precision matrix, which means that some edges will be zero. A natural way 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 penalized log-likelihood that weighs the variables by their scale: over positive definite symmetric matrices, where τ ≥ 0 is a penalty parameter for the offdiagonal elements. We refer to the objective function in (2.1) as a "weighted penalized log-likelihood". The subscript τ in Θ τ means that the solution of the optimization problem in (2.1) will depend upon the choice of the tuning parameter. In order to simplify notation, we will omit the subscript.
One of the popular and fast algorithms to solve the optimization problem in (2.1) is called the Graphical LASSO (GLASSO), which was introduced by Friedman et al. (2008). Define the following partitions of W, S and Θ: Let β ≡ −θ 12 /θ 22 . The idea of GLASSO is to set W = S + τ I in (2.1) and combine the gradient of (2.1) with the formula for partitioned inverses to obtain the following 1regularized quadratic program β = arg min over positive definite symmetric matrices. 3: Repeat for j = 1, . . . , p, 1, . . . , p, . . . until convergence: • Partition W into part 1: all but the j-th row and column, and part 2: the j-th row and column.
• Solve the score equations using the cyclical coordinate descent: This gives a (p − 1) × 1 vector solution β.
• Update w 12 = W 11 β. As was shown in Friedman et al. (2008), the estimator produced by Algorithm 1 is guaranteed to be positive definite. Furthermore, Janková and van de Geer (2018) showed that Algorithm 1 is guaranteed to converge and produces consistent estimator of precision matrix under certain sparsity conditions.

Approximate Factor Models for Forecast Errors
The approximate factor models for the forecasts were first considered by Chan et al. (1999). 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 1 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 Diebold and Shin (2019), 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 1, 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, . . . , T and e t = (e 1t , . . . , e pt ) ∼ N (0, Σ) is a p × 1 vector of forecast errors. Assume that the generating process for the forecast errors follows a q-factor model: . . , f qt ) are the common factors of the forecast errors for p models, B is a p × q matrix of factor loadings, and ε t is the idiosyncratic component that cannot be and Watson (2002). Strict factor structure assumes that the idiosyncratic forecast error terms, ε t , are uncorrelated with each other, whereas approximate factor structure allows correlation of the idiosyncratic components (Chamberlain and Rothschild (1983)).
We use the following notations: be the precision matrices of forecast errors, idiosyncratic and common components respectively. The objective function to recover factors and loadings from (3.1) is: It is well-known (see Stock and Watson (2002) among others) that B estimated from the first q eigenvectors of 1 T T t=1 e t e t is the solution to (3.4). Given a sample of the estimated Moving forward to the forecast combination exercise, suppose we have p competing forecasts, y t = (ŷ 1,t , . . . ,ŷ p,t ) , of the variable y t , t = 1, . . . , T . The forecast combination is defined as follows: where w is a p × 1 vector of weights. Define a measure of risk MSFE(w, Σ) = w Σw. As shown in Bates and Granger (1969), the optimal forecast combination minimizes the MSFE of the combined forecast error: where ι p is a p × 1 vector of ones. The solution to (3.6) yields a p × 1 vector of the optimal forecast combination weights: If the true precision matrix is known, the equation (3.7) guarantees to yield the optimal forecast combination. In reality, one has to estimate Θ. Hence, the out-of-sample performance of the combined forecast is affected by the estimation error. As pointed out by Smith and Wallis (2009), when the estimation uncertainty of the weights is taken into account, there is no guarantee that the "optimal" forecast combination will be better than the equal weights or even improve the individual forecasts. Define a = ι p Θι p /p and a = ι p Θι p /p. We can 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 Θ. More details are discussed in section 4 and Theorem 1.

Factor Graphical LASSO for Forecast Errors
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 Koike (2020), 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. Hence, instead of imposing sparsity assumption on the precision of forecast errors, Θ, we require sparsity of the precision matrix of the idiosyncratic errors, Θ ε . The latter is obtained using the estimated residuals after removing the comovements induced by the factors (see Brownlees et al. (2018); Koike (2020)). Naturally, once we condition on the common components, it is sensible to assume that many remaining partial correlations of ε t will be negligible and thus Θ ε is sparse.
We use the weighted Graphical LASSO as a shrinkage technique 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 To get Θ ε , we use the weighted GLASSO Algorithm 1, with the initial estimate of the covariance matrix of the idiosyncratic errors Once we estimate Θ f and Θ ε , we can get Θ using a sample analogue of (4.1). We call the proposed procedure Factor Graphical LASSO and summarize it in Algorithm 2.
Step 1 and Θ ε from Step 2 to estimate Θ using the sample counterpart of the Sherman-Morrison-Woodbury formula in (4.1): Let Θ ε,τ be the solution to (4.2) for a fixed τ . To choose the optimal shrinkage intensity coefficient, we minimize the following Bayesian Information Criterion (BIC) using grid search: The grid G ≡ {τ 1 , . . . , τ M } is constructed as follows: the maximum value in the grid, τ M , is set to be the smallest value for which all the off-diagonal entries of Θ ε,τ M are zero, that is, the maximum modulus of the off-diagonal entries of Σ ε . The smallest value of the grid, τ 1 ∈ G, is determined as τ 1 ≡ ϑτ M for a constant 0 < ϑ < 1. The remaining grid values τ 1 , . . . , τ M are constructed in the ascending order from τ 1 to τ M on the log scale: We use ϑ = log p/T + 1/ √ p (motivated by the convergence rate from Theorem 1) and M = 10 in the simulations and the empirical exercise.
We can use Θ to estimate the forecast combination weights w where Θ is obtained from Algorithm 2. 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.

Asymptotic Analysis of the Estimator
We first introduce some terminology and notations. Let A ∈ S p . Define the following set for j = 1, . . . , p: And assume that (ii) ι p Θι p /p ≥ c > 0, where c > 0 is a positive constant.
(A.2) (Pervasive factors) There exists a positive definite q × q matrixB such that We also impose strong mixing condition.
Assumption (A.1) 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 (A.2). 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 (A.1) is needed for consistent estimation of the optimal forecast combination weights. Assumptions Let Λ q = diag(λ 1 , . . . , λ q ) be a matrix of q leading eigenvalues of Σ, and V q = (v 1 , . . . , v q ) is a p × q matrix of their corresponding leading eigenvectors. Define Σ, Λ q , V q to be the estimators of Σ, Λ q , V q . We further let Λ q = diag(λ 1 , . . . ,λ q ) and V q = (v 1 , . . . ,v q ) to be constructed by the first q leading empirical eigenvalues and the corresponding eigenvectors Fan et al. (2018), we require the following bounds on the componentwise maximums of the estimators: Assumptions ( In addition, the following structural assumptions on the model are imposed: Note that Assumptions (B.1)-(B.3) and (C.1) are standard assumptions and are used in Fan et al. (2018).
To study the properties of the combination weights in (4.5) and MSFE, we first need to establish the convergence properties of precision matrix produced by Algorithm 2. Let . The deterministic sequences s T and d T will control the sparsity Θ ε for Factor GLASSO. Note that d T can be smaller than or equal to s T .
Let 1T be a sequence of positive-valued random variables such that −1 1T ω T P − → 0 and 1T d T s T P − → 0, with τ ω T (where τ is the tuning parameter for the Factor GLASSO in (4.2)). Lee and Seregina (2020)  and Factor GLASSO.
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 1 can be found in Appendix A. Note that the rates of convergence for MSFE and precision matrix Θ are the same and both are faster than the combination weight rates. In contrast to the classical graphical model in Algorithm 1, the convergence properties of which were examined by Janková and van de Geer (2018) among others, the rates in Theorem 1 depend on the sparsity of Θ ε rather than of Θ. This means that instead of assuming that many partial correlations of forecast errors e t are negligible, which is not realistic under the factor structure, we impose a milder restriction requiring many partial correlations of ε t to be negligible once the common components have been taken into account.

RD-Factor GLASSO for Forecast Errors
There are two streams of literature that study time-varying networks. The first one models dynamics in the precision matrix locally. Zhou et al. (2010)  We augment the framework in Section 4 to account for regime switching by modeling the change in precision matrix due to N known structural breaks. Define n i ≡ t i − t i−1 to be the sample between the i-th and (i − 1)-th break points, where i = 1, . . . , N , N i=1 n i = T , We propose that the dynamics of the system evolves through the precision matrix of the idiosyncratic component. Let Σ ε,i and Σ i be covariance matrices of idiosyncratic part and forecast errors in regime i. Define the corresponding precision matrices to be Θ ε,i ≡ Σ −1 ε,i and We assume Σ f i = Σ f for all regimes i for simplicity but it can be generalized.
we use the following opti-mization problem: where the penalty for the off-diagonal (od) elements is Θ ε,i od,1 = l =q γ ε,ll,i γ ε,qq,i |θ ε,lq,i |, γ ε,ll,i is the (l, l)-th element of Γ 2 ε,i ≡ diag( Σ ε,i ) and θ ε,lq,i is the lq-th element of matrix Θ ε,i . Figure 2 visualizes dynamics of the precision matrix. Figure 2: Change of precision matrix over time: β is the penalty that enforces temporal consistency and ψ is a convex penalty function.
The optimization problem in (5.1) has two tuning parameters: α, which determines the sparsity level of the network, and β, which controls the strength of resemblance between two neighboring precision estimators. In simulations and the empirical application we use the following procedure for tuning α and β: first, we set a grid of values (α, β) ∈ {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 α and β in the remaining 1/3 to yield the smallest value of the objective function, which is chosen to be either |||·||| 2 -loss of precision matrix for simulations in Subsection 6.1, or MSFE for simulations in Subsection 6.2 and the empirical application. Note that when β = 0, the optimization in (5.1) reduces to estimating Θ ε,i using Algorithm 2 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 β is large, there are weak structural breaks in Θ ε,i , and Θ ε,i 's are estimated by using the data across different regimes. Section 7 provides more discussion on this in the context of our empirical application.
The smoothing function ψ(·) in (5.1) can be LASSO (ψ = l,q |·|), Group LASSO (ψ = q · q 2 ), Ridge (ψ = l,q (· lq ) 2 ), Max norm penalty (ψ = q max l |· lq |). 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. For max-norm penalty, if an element of the precision matrix changes by in regime i, then other elements are allowed to change by up to without any additional penalty. This penalty allows to capture changes in the clusters of forecast errors, whereas keeping the rest of the network unchanged. In our empirical application we use Ridge penalty to accommodate smooth transitions of precision over time.
To estimate (5.1) we use the alternating direction method of multipliers (ADMM) algorithm described in details in Supplementary Appendix B. Once Θ ε,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 Regime Dependent Factor GLASSO (RD-Factor GLASSO) and summarize it in Algorithm 3.

Algorithm 3 RD-Factor GLASSO
1: Estimate f t and B t in (3.1). Get Σ f , Θ f and ε t = B f t − e t .
3: Use Θ ε,i , Θ f and B from Steps 1-2 to get Θ i : We develop a scalable implementation of (5.1) for the RD-Factor GLASSO in Algorithm 3 through ADMM, which is extensively discussed in Supplementary Appendix B.

Monte Carlo
We divide the simulation results into two subsections. In the first subsection we study the consistency of the Factor GLASSO and RD-Factor GLASSO 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 2, 3 with equal-weighted (EW) forecast combination, and combinations that use GLASSO without factor structure (Algorithm 1). Similarly to the literature on graphical models, all exercises use 100 Monte Carlo simulations.
We show consistency of the Factor GLASSO in (i) the operator norm, Θ − Θ 0 2 , and (ii) in 1 -vector norm for the combination weights, w − w 1 , where w is given by (3.7).
The forecast errors are assumed to have the following structure: where e t is a p × 1 vector of forecast errors following N (0, Σ), f t is a q × 1 vector of factors, B is a p × q matrix of factor loadings, φ f is an autoregressive parameter in the factors which is a scalar for simplicity, ζ t is a q × 1 random vector with each component independently following N (0, σ 2 ζ ), ε t is a p × 1 random vector following N (0, Σ ε ), with sparse Θ ε that has a random graph structure described below. To create B in (6.1) we take the first q columns of an upper triangular matrix from a Cholesky decomposition of the p × p Toeplitz matrix parameterized by ρ: that is, . . , p}. We set ρ = 0.2, φ f = 0.2 and σ 2 ζ = 1. The specification in (6.1) leads to the low-rank plus sparse decomposition of the covariance matrix: When Σ ε has a sparse inverse Θ ε , it leads to the low-rank plus sparse decomposition of the precision matrix Θ, such that Θ can be expressed as a function of the low-rank Θ f plus sparse Θ ε .
We consider the following setup: let p = T δ , δ = 0.85, q = 2(log(T )) 0.5 and T = [2 κ ], for κ = 7, 7.5, 8, . . . , 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 Θ ε is constructed as follows: we first generate the adjacency matrix using a random graph structure. Define a p × p adjacency matrix A ε which represents the structure of the graph: where a ε,ij denotes the i, j-th element of the adjacency matrix A ε . We set a ε,ij = a ε,ji = 1, for i = j with probability π, and 0 otherwise. Such structure results in s T = p(p − 1)π/2 edges in the graph. To control sparsity, we set π = 500/(pT 0.8 ), which makes s T = 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 (Pedregosa et al. (2011)). 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. Figure 3 shows the averaged (over Monte Carlo simulations) errors of the estimators of the precision matrix Θ 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 Ledoit and Wolf (2004)). As evidenced by Figure 3, Factor GLASSO demonstrates superior performance over EW and non-factor based model (GLASSO). Furthermore, our method achieves lower estimation error in the combination weights (3.9), which leads to lower risk of the combined forecast as shown in (3.8). 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 good convergence properties. To incorporate structural break in Θ ε we add the following modification to the DGP described above. 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 Θ ε,1 , the value of the largest coefficient is set to 0.4; whereas in the precision matrix of the idiosyncratic errors after the break, Θ ε,2 , the value of the largest coefficient is set to 0.6. As a consequence, even though both matrices are still sparse, Θ ε,2 has larger partial correlations. We use Θ ε,1 and Θ ε,2 to generate ε t in (6.1). Figure 4 shows the performance of all models including RD-Factor GLASSO: accounting for the break in Θ ε significantly reduces the estimation error of precision matrix and combination weights. Notice that Factor GLASSO still outperforms GLASSO, although the dominance is less pronounced compared to the setup without structural breaks.

Comparing Performance of Forecast Combinations
We consider the standard forecasting model in the literature (e.g., Stock and Watson (2002)), 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, x t is an N × 1 vector of regressors (predictors), β is an N × 1 parameter vector, g t is an r × 1 vector of factors, Λ is an N × r matrix of factor loadings, v t is an N × 1 random vector following N (0, σ 2 v I N ), φ is an autoregressive parameter in the factors which is a scalar for simplicity, ξ t is an r × 1 random vector with each component independently following N (0, σ 2 ξ ), t+1 is a random error following N (0, σ 2 ), and α is an r × 1 parameter vector which is drawn randomly from N (1, 1). We set σ = 1. The coefficients θ s are set according to the rule θ s = (1 + s) c 1 c s 2 , (6.8) as in Hansen (2008). We set c 1 ∈ {0, 0.75} and c 2 = 0.9. We set N = 100 and generate r = 5 factors. To create Λ in (6.5) we take the first r rows of an upper triangular matrix from a Cholesky decomposition of the N ×N Toeplitz matrix parameterized by ρ = 0.9. The ranking of competing models was not very sensitive to varying values of φ, ρ, c 2 , and r -the results examining sensitivity to a grid of 10 different AR(1) coefficients φ equidistant between 0 and 0.9, a grid of 10 different values of ρ equidistant between 0 and 0.9, c 2 ∈ {0.6, 0.7, 0.8, 0.9}, and r = [1, · · · , 7] are available upon request.
The total number of observations is T . The period for training the models is set to be m 1 = T /2 -this is used to train competing FAR models in (6.9). The remaining part of the sample, m 2 = T /2 is split as follows: the estimation window for training competing models (that is, EW, GLASSO, Factor GLASSO, and RD-Factor GLASSO) is set to be window = 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 t = 1, . . . , m/2. Recall that q denotes the number of factors in the forecast errors as in equation (3.1). Figure 5 shows the MSFE for different sample sizes and fixed parameters: we report the results for two values of c 1 ∈ {0, 0.75}. As evidenced from Figure 5, the models that use the factor structure outperform EW combination and non-factor based counterparts for both values of c 1 .
To incorporate structural break we add the following modification to the DGP described above. The period for training the models is set to be m 1 = T /3 -this is used to train competing FAR models in (6.9). The remaining part of the sample, m 2 = 2 × T /3 is split as follows: the estimation window for training competing models is set to be window = m 2 /2. Figure 5: Plots of the MSFE over the sample size T . c 1 = 0 (left), c 1 = 0.75 (right), c 2 = 0.9, N = 100, r = 5, σ ξ = 1, L = 7, K = 2, p = 24, q = 3, ρ = 0.9, φ = 0.8.
We roll the estimation window over the 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 θ s in (6.8) we set c 2 = 0.3, and after the break c 2 = 0.9. All other parameters stay unchanged. Figure 6 shows the performance of all models including RD-Factor GLASSO: similarly to the conclusions in the previous subsection, accounting for the break significantly reduces MSFE of the combined forecast.

Application to Combining ECB SPF Forecasts
We use quarterly forecasts on the expected rates of inflation, real GDP growth and unemployment rate in the Euro area published by the ECB. The raw data records 119 Figure 6: Plots of the MSFE over the sample size T . c 1 = 0.75, c 2 = 0.3 (before the break), c 2 = 0.9 (after the break), N = 100, r = 5, σ ξ = 1, L = 7, K = 2, p = 24, q = 3, ρ = 0.9, φ = 0.8. 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 et al. (2020) to obtain most qualified forecasters: first, we filter out irregular respondents if they missed more than 50% of the observations; second, we use a random forest imputation algorithm (Stekhoven (2022); Stekhoven and Buhlmann (2012) We consider three choices of the training sample: m ∈ {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 (3.1)) is chosen using the standard data-driven method that uses the information criterion IC1 described in Bai and Ng (2002). In the majority of the cases the optimal number of factors was estimated to be equal to 1. To explore the benefits of using Factor GLASSO and RD-Factor GLASSO for forecast error quantification, we consider several alternative estimators of covariance/precision matrix of the idiosyncratic component in (4.1): (1) linear shrinkage estimator of covariance developed by Ledoit and Wolf (2004) further referred to as Factor LW; (2) nonlinear shrinkage estimator of covariance by Ledoit and Wolf (2017) (Factor NLW); (3) POET (Fan et al. (2013)); (4) constrained 1 -minimization for inverse matrix estimator, CLIME (Cai et al. (2011)) (Factor CLIME); (5) nodewise regression developed by Meinshausen and Bühlmann (2006) (Factor MB). To examine the benefits of imposing sparsity on Θ ε we also include the factor model without sparsity assumption on the idiosyncratic error precision matrix (referred to as Not Sparse) -this corresponds to imposing τ = 0 in (4.2). Our benchmark is the simple average with equal weights on all forecasters (referred to as EW). Going back to the discussion in Section 5 regarding setting β = 0 in equation (5.1): as we pointed out, this corresponds to using only post-break sample for estimation which is suboptimal since the value of β is already chosen optimally from the grid that includes β = 0 to minimize the MSFE. Hence, by construction, RD-Factor GLASSO is superior to using only post-break data.
For RD-Factor GLASSO we consider the following structural break points associated with a global financial crisis of 2007-09: for real GDP series the break is taken to be 2009:Q1, for inflation series the break is 2008:Q2, and for unemployment rate -2008:Q2. We chose these points based on the behavior of the actual series which exhibited a change in volatility following the break. This corresponds to one structural break for each series, or, in terms of the notations used in (5.1), N = 2 (two regimes). Table 1 compares the performance of Factor GLASSO and RD-Factor GLASSO with the competitors for predicting three macroeconomic indicators for Euro-area using a combination of ECB SPF forecasts. For the inflation series training the models is computationally challenging due to p ≥ T , hence the numerical routine used in several models (Factor NLW and Factor CLIME) failed to converge and we were not able to get the final estimates despite making appropriate changes to the tuning parameters used in NLW and CLIME. Furthermore, GLASSO also had difficulty in handling inflation series since for all m almost all off-diagonal elements of Θ were shrunk to zero. As a consequence, for the inflation series GLASSO produced the same results as EW.
There are three main findings that we learn from analyzing Table 1: (1) for all series factor-based models outperform non-factor ones (such as EW and GLASSO). This means that incorporating the factor structure in the forecast errors improves forecasting performance.
(2) for all series the Not Sparse model provides the worst performance. 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) RD-Factor GLASSO is the best-performing model for real GDP and inflation series, and Factor GLASSO is the second-best performing model. For the unemployment rate, Factor GLASSO outperforms RD-Factor GLASSO. This result is supported by the behavior observed in the actual series: real GDP and inflation exhibit a clear break in variance following the global financial crisis of 2007-09, however this is not the case for the unemployment rate series that did not have noticeable changes in variance throughout the whole sample period.

Conclusions
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.
Our first algorithm, Factor GLASSO, is shown to consistently estimate forecast combination weights and MSFE. For the ease of practical use we develop a scalable optimization procedure for RD-Factor GLASSO, 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.
where the last rate is obtained using the assumptions of Theorem 1.

A.1 Proof of Theorem 1
First, note that the forecast combination weight can be written as As shown in Callot et al. (2019), the above can be rewritten as Prior to bounding the terms in (A.3), we first present an inequality which is used in the derivations. Let A ∈ R p×p and v ∈ R p×1 . Also, let A j and A j be a p × 1 and 1 × p row and column vectors in A, respectively.
Hölders inequality was used to obtain each inequality in (A.4). If A ∈ S p , then the last expression can be further reduced to p|||A||| 1 v ∞ . Let us now bound the right-hand side of (A.3). In the numerator we have: the rates was derived in Lee and Seregina (2020), and the inequality follows from (A.4).
where the rate follows from Lemma 1 (a) and the inequality is obtained from (A.4). Combining (A.5), (A.6), and Lemma 1 (c) we get: where the last equality holds under the assumptions of Theorem 1. For the denominator of (A.3) it easy to see that | a|a = O P (1) using the results of Lemma 1 (b).
For the MSFE part of Theorem 1, using Lemma 1 (b)-(c), we get where the last rate is obtained using the assumptions of Theorem 1.

Appendix B Implementation via ADMM Algorithm
To enable practical implementation of the RD-Factor GLASSO, we develop an optimization procedure using ADMM algorithm to solve the convex optimization problem in (5.1).
We will solve a separate update for each (Z i,2 , Z i−1,1 ) pair for i = 2, . . . , N : (Z k+1 i,2 , Z k+1 i−1,1 ) = arg min Note that (B.9) is guaranteed to converge to a fixed point since it can be written as a proximal operator: (Z k+1 i,2 , Z k+1 i−1,1 ) = prox β ρ ψ(·) Θ ε,i + U i,2 , Θ ε,i−1 + U i−1,1 (B.10) Remark 1. A proximal operator of the scaled function νf , where ν > 0 can be expressed as: where f is a closed proper convex function. Note: Parikh and Boyd (2014) show that the fixed points of the proximal operator of f are precisely the minimizers of f , i.e. prox νf (x ) = x if and only if x minimizes f .

The Θ step:
The updating rule in (B.5) can be further simplified to obtain a closed-form solution. Rewrite (B.5): , and η = n i 3ρ .
Take the gradient of the updating rule in (B.11) in order to get an analytical solution: Equation (B.13) implies that Θ k+1 ε,i and 1 η A k − Σ ε,i share the same eigenvectors.
Let Q i Λ i Q i be the eigendecomposition of 1 η A k − Σ ε,i , where Λ i = diag(λ 1,i , . . . , λ p,i ), and Q i Q i = Q i Q i = I. 1 Pre-multiply (B.13) by Q i and post-multiply it by Q i : Now construct a diagonal solution of (B.14): (B.15) whereṽ i,j denotes the j-th eigenvalue of Θ ε,i . Solving forṽ j,i we get: v j,i = λ j,i + λ 2 j,i + 4 η 2η −1 . (B.16) Now we can calculate Θ k+1 ε,i which satisfies the optimality condition in (B.14): Use the definition of η = n i 3ρ : Step (B.18) is the most computationally intensive task in the algorithm since the runtime of decomposing a p × p matrix is O(p 3 ). Also, note that compared to standard ADMM without smoothing penalty β, (B.18) enforces stronger shrinkage. This is consistent with our motivation for the additional constraint -to smooth the estimator of precision matrix.