Bootstrap VAR forecasts: The effect of model uncertainties

VAR models are popular to forecast macroeconomic time series. However, the model, the parameters, and the error distribution are rarely known without uncertainty, so bootstrap methods are applied to deal with these sources of uncertainties. In this paper, the performance of the popular forecast Bonferroni cubes based on the Gaussian method and variants of the bootstrap procedure that incorporate error distribution, parameter uncertainty, bias correction, and lag order uncertainty are compared. Monte Carlo simulations suggest that the best performance of bootstrap cubes are obtained when the parameter uncertainty is considered, being the bias and model uncertainties important for long-run forecast regions in persistent VAR models. Similar conclusions are found in an empirical application based on a four variate system containing US monthly percent changes of the industrial production index, the S&P500 stock market index, its dividend yield, and the unemployment rate.


| INTRODUCTION
When forecasting the evolution of multivariate time series systems, analysts are increasingly interested in obtaining the full probability distribution; see, for example, Diebold et al. (1999), Jore et al. (2010) and Giacomini and Ragusa (2014).A popular model to forecast macroeconomic and financial time series is the VAR(p) model; see, for example, Stock and Watson (2001).The standard procedure to construct VAR forecast densities assumes Gaussian errors and known lag order and parameters; see, for instance, Lütkepohl (1991).However, in practice, the lag order needs to be selected, the parameters need to be estimated, and the error distribution is rarely known.Consequently, the forecast model is actually an estimated VAR model that is treated as if it were the true data generating process (DGP).This is to say that forecast densities are conditional on the selected lag order, the estimated parameters, and the errors distribution with uncertainties due to these factors being ignored.Furthermore, the least squares (LS) parameter estimator can suffer from small-sample bias for processes that have roots close to the unit circle so density forecasts may be contaminated; see Lawford and Stamatogiannis (2009) for non-stationary VAR models.
Some of these issues can be addressed by the asymptotic theory.This is the case of the bias of the LS estimator, which can be adjusted using, for instance, the correction proposed by Pope (1990).Likewise, the parameter uncertainty can be incorporated into the forecast variability by using the asymptotic distribution of the LS estimator as proposed by Reinsel (1980) and Lütkepohl (1991).The lag order uncertainty can be selected by an information criteria, but as far as we know, there is no asymptotic proposal to incorporate the uncertainty surrounding this selection into forecast densities; see Schorfheide (2005) for an analysis of the effect of the lag order uncertainty on the forecast uncertainty.Alternatively, Bayesian procedures can generate forecast densities that incorporate the model and parameter uncertainties; see Koop (2013) and Cross et al. (2020) for some references.Notwithstanding, Bayesian methods are strongly tied to the error Gaussianity assumption because their computation burden is very strong when this assumption is abandoned.
Bootstrap procedures can deal with all the sources of model uncertainty described above.Thus, several authors have pointed out the attractiveness of bootstrap methods to forecast VAR models because they are able to incorporate lag order and parameter uncertainties without assuming any particular error distribution; see Berkowitz and Kilian (2000) and Cavaliere et al. (2020) for surveys.This is the reason why bootstrap forecast densities and their associated forecast regions have often been considered for VAR models; see, for example, Staszewska-Bystrova and Winker (2013), Wolf and Wunderli (2015), Fresoli et al. (2015) and Beyaztas (2019) for some applications.However, there have been no attempts to simultaneously consider the different sources of uncertainty to assess their relative importance.
By building bootstrap regions, this paper aims to analyze the contribution of the model parameter, its bias, the lag order, and error uncertainties in the context of multivariate VAR(p) models.To measure the contribution of each source of uncertainty, we generate artificial time series of different sizes from VAR models with different number of parameters, autoregressive roots, and error distributions.First, we construct forecast regions using the standard approach based on estimated VAR models with Gaussian errors.Second, we implement the bootstrap procedure proposed by Fresoli et al. (2015) incorporating the uncertainty about the error distribution with fixed parameters and lag order.Third, we consider a bootstrap procedure that deals with the parameter uncertainty.Then, we add bias-correction.Finally, we implement the bootstrap incorporating also the lag order uncertainty.We compare the alternative procedures in terms of the empirical coverages of the Bonferorri forecast cubes which are often implemented when constructing forecast regions; for instance, see Jord a et al. (2013)., Staszewska-Bystrova and Winker (2014) and Beyaztas (2019).The conclusions of the Monte Carlo experiments are important for practitioners given that the relative importance of each source of uncertainty is established in terms of the model dimension, error distribution, persistence of shocks, and sample size.To illustrate the relevance of the results for real systems of time series, the proposed procedures are implemented to obtain forecast densities of a macroeconomic system.
There are procedures to construct confidence regions with some nominal coverage for a multivariate system at a given forecast horizon (i.e., forecast regions).For instance, we may construct ellipsoids with a given nominal coverage, but they are based on a Gaussianity assumption and are hard to obtain for high dimensional systems; see, for instance, Kim (1999).Alternatively, we can obtain Bonferroni cubes with a lower bound nominal coverage that are easy to implement as they are formed by using individual forecast intervals for each of the variables in the system; see Fresoli et al. (2015) and Wolf and Wunderli (2015).In this paper, instead of pursuing a comparison across different forecast regions, we aim to assess how different sources of uncertainty affect the empirical coverages of the Bonferroni forecast cubes.It is also worth mentioning the growing literature that analyzes methods to construct confidence paths, either for impulse response coefficients or out-of-sample forecasts; see, for example, Lütkepohl et al. (2015), Bruder and Wolf (2018) and Lütkepohl et al. (2020).Identifying and assessing the different sources of uncertainty are critical to the construction of forecast regions or paths with a certain nominal coverage.
The structure of this paper is the following.Section 2 establishes notation.By describing briefly the VAR(p) model and the construction of forecast densities and regions using the standard Gaussian and bootstrap procedures, we provide in Section 2 all the necessary notation.Section 3 describes the Monte Carlo study based on VAR models with different orders, persistence properties and error distributions.Section 4 contains an empirical application to a four-dimensional macroeconomic system.Finally, Section 5 presents the conclusions.

| FORECASTING WITH VAR MODELS
Consider the following multivariate VAR(p) model where y t is the N Â 1 vector of observations at time t, c is a N Â 1 vector of constants, Φ j are N Â N autoregressive matrices, for j ¼ 1, …, p, satisfying the stationary conditions and ε t is a sequence of independent white noise N Â 1 vectors with nonsingular contemporaneous covariance matrix Σ ε .Let Π ðpÞ ¼ Φ 1 , …,Φ p Â Ã be the N Â Np matrix of autoregressive parameters.Stationarity guarantees that the VAR process has a vector moving average (VMA) representation given by y Denote by ĉ and ΠðpÞ the matrices containing the LS estimates of the parameters obtained after selecting the lag order p by the Akaike information critetia (AIC) or Schwartz criteria (BIC), and by Σε and Ψi the sample covariance matrix of the residuals and the estimates of the MA matrices, respectively.Then, the standard hsteps-ahead forecast density is constructed-after assuming Gaussian errors-as follows where ŷTþhjT ¼ ĉ þ P p i¼1 Φj ŷTþhÀijT , with ŷTþjjT ¼ y Tþj for j < 0, is the conditional mean of y T + h given {y 1 , … , y T }, which is obtained by substituting the unknown lag order and parameters by their corresponding estimates, and Σy ðhÞ ¼ P hÀ1 i¼1 Ψi Σε Ψ0 i .The forecast densities in (2) are denoted as standard (STD).The STD densities can be inadequate to construct forecast regions for several reasons.First, they do not tackle the parameter bias of the LS estimates which affects both ŷTþhjT and Σy ðhÞ.Second, Σy ðhÞ partially reflects the uncertainty around the point forecast because other sources of uncertainty, such as those associated with the selection of the lag order and the estimation of the parameters, are omitted.Finally, the Gaussianity assumption may be misleading for many macroeconomic and financial systems; see Kilian (1998b).The bias of the LS estimator can be corrected by using an asymptotic bias-correction formula as proposed by Pope (1990).The matrix of the corrected LS estimates Πc ðpÞ will be used to obtain ŷc TþhjT and Σc y ðhÞ.Bootstrap procedures implemented in this paper to construct multivariate forecast densities are based on the proposal in Fresoli et al. (2015).This procedure is of interest because it does not use the backward representation (BR) of the VAR models as needed by, for example, Kim (1999); Kim (2004), Grigoletto (2005) and Staszewska-Bystrova (2011).Consequently, it is simple and widely applicable, making it also attractive to forecast more complex models; see Fresoli and Ruiz (2016) and Bauwens et al. (2017) for similar procedures to forecast dynamic conditional correlations and realized covariance matrices, respectively.
Steps involved in the bootstrap procedure to construct h-steps-ahead forecast densities for VAR(p) models are described in detail in Fresoli et al. (2015).By applying this method, we obtain fŷ * ð1Þ TþhjT , …, ŷ * ðBÞ TþhjT g, with B being the number of bootstrap replicates, which can be used to approximate the h-steps-ahead forecast distribution of the process.More specifically, if G * ðxÞ ¼ Pðy * Tþh < xÞ is the distribution function of y * Tþh , then its Monte Carlo estimate is given by Ĝ where I(Á) is an indicator function.The estimated bootstrap distribution can be used either to construct forecast intervals and regions with appropriate probability content or to provide the probability forecasts associated to some events of interest; see Garratt et al. (2003).The asymptotic validity of the procedure is analyzed by Fresoli et al. (2015).
When the bootstrap procedure is run with a fixed lag order and given parameters, using the parameter estimates without bias-correction, the bootstrap forecast densities (denoted as distribution bootstrap, DB) only incorporate errors distribution uncertainty; for example, STATA implements a somewhat similar approach to DB in which out-of-sample bootstrap draws are constructed conditioned on estimated parameters, lag order, and the last p observations (StataCorp, 2013). 1 Second, we consider the bootstrap procedure that avoids bias-correction and does not re-estimate the lag order for each of the bootstrap series.This bootstrap procedure, which is called parameter bootstrap (PB), tackles, in addition to the error distribution uncertainty, the forecast variability due to parameter uncertainty; see, for instance, Wolf and Wunderli (2015) who construct bootstrap forecast regions using this procedure.Third, the PB procedure can be implemented with bias-corrected parameters; see Clements and Kim (2007) who highlight that in the presence of biased parameters, the resulting bootstrap forecast regions may be double biased because the bootstrap estimates of the parameters might also suffer from bias.This procedure implements the asymptoticbias correction formula of Pope (1990) and is denoted as bias-corrected parameter bootstrap(BCB). 2A bootstrap procedure of this type is the one proposed by Alonso et al. (2018) for forecasting highly persistent dynamic autoregressive factors.Finally, we consider the 1 The difference with respect to DB is that, however, STATA forecast intervals are obtained using the standard formula ŷ * i, Tþh AE z 1À α 2 Σ * y i ðhÞ with the ith individual point and MSFE replaced by their bootstrap counterparts.

2
There are alternatives to the asymptotic bias-correction formula as the bootstrap-after-bootstrap or doubled bootstrap approach; see Kilian (1998c).However, the latter does not guarantee a better performance than the asymptotic formula in terms of empirical coverages but is computationally more demanding.bias-corrected parameter and endogenous bootstrap (EB) in which the uncertainty about the lag order is tackled by re-estimating it for each bootstrap series; see, for example, Staszewska-Bystrova and Winker (2013) and Schreiber and Soldatenkova (2016).
It is worth noting that finite order VAR can be considered as an approximation to finite order VMA, infinite order VAR, or VAR models exhibiting some non-linearity (e.g., multivariate GARCH errors).Model uncertainty arising from the approximation of a general DGP with a finite order VAR model is not addressed in this paper.

| MONTE CARLO EXPERIMENTS
To asses the impact of bias and model uncertainty on multivariate Bonferroni bootstrap forecast regions, we carry out Monte Carlo experiments in this section.We consider six bivariate VAR(p) models with lag orders p ¼ 2, 4, and 8, and two VAR(2) models with four variates.For each model, we choose autoregressive matrices with different persistence dynamics. 3Table 1 summarizes the eight models considered.Four error distributions are assumed, namely, the Gaussian, Student-3, χ 4 , and Bimodal 0.7ÂN(1, 1) + 0.3 Â N(À 7/3, 1) distributions.We generate multisame covariance matrix, given by vechðΣ ε Þ ¼ ð1:84, 0:03, 0:03Þ and (1.84, 0.03, 0.03, 0.00, 0.03, 0.00, 0.00, 0.02, 0.002, 0.03) for the two and four variable models, respectively, where vech denotes the lower diagonal column stacking operator.M ¼ 2000 replicates of sizes T ¼ 100 and 250 are generated for each of the models considered.The sample sizes assumed here are representative of those commonly used in practice and our design includes models with different VAR dimension, model lag, autoregressive roots, and error distributions so that their effects can be explored.variateerrors by drawing samples from these independent univariate distributions.Errors are then centered and appropriately re-scaled through a Cholesky matrix to have the For each replicate, AIC or BIC is implemented to choose the lag order, with the maximum order being equal to 8, 12, and 16 for the models with p ¼ 2, 4, and 8, respectively.After estimating the parameters by LS, h-steps-ahead forecast densities for h ¼ 1, …, 12 are constructed assuming Gaussian errors (STD) as in (2).Furthermore, bootstrap forecast densities are constructed based on B ¼ 2000 bootstrap replicates by implementing the DB, PB, BCB and EB procedures.In each case, we construct the corresponding 100Â(1-α)% Bonferroni cube as follows: where q i (τ) is the τth quantile of the empirical marginal forecast distribution of ŷ * i,Tþh .We consider 100α%=5% and 10%.After generating 2000 future true values of y T + h , and counting the number of them lying inside the 90% and 95% forecast Bonferroni cubes, we obtain the corresponding empirical coverages that are used to asses the adequacy of the STD, DB, PB, BCB, and EB densities.
Figure 1 reports the differences between the empirical and nominal coverages of the Bonferroni cubes for h ¼ 1 and 12, for the different procedures and VAR models when T ¼ 100, the nominal coverage is 90% and the lag order is estimated by AIC.The first row displays results when the errors are Gaussian.Note that, in this case, the STD cubes are constructed using the true forecast error distribution.For h ¼ 1, STD have coverages that are always below the nominal.The undercoverage can be severe, reaching, for instance, 8% for a highly persistent model like M6 in which p ¼ 8.Note that, in models with smaller persistence, the coverages are closer to the 3 In Appendix 2, we report the autoregressive matrices for each of the eight VAR models used in this MC experiments.nominal when h ¼ 12, but in those with higher persistence, the coverages are worse as the forecast horizon increases.If we look at the empirical coverages of DB, it is observed that, for h ¼ 1 and T ¼ 100, they are slightly smaller than those of STD regions, regardless of the VAR models under consideration.Although this result could be expected given that the true error distribution is Gaussian, it is important to point out that for h ¼ 12 the coverages of the DB cubes are comparable to those of the STD cubes.By including the parameter uncertainty using PB, we obtain Bonferroni regions with accuracy closer to the nominal, for all VAR models if T ¼ 100 and h ¼ 1.The improvement is larger for large and high persistent VAR models as, for example, model M6, and h ¼ 12; see Kim (1999) for a similar result.
With respect to the VAR dimension, note the differences between the empirical and nominal coverages of the STD cubes as we move from M1 to M8.For instance, increasing the number of variables to four-as in models M7 and M8-produces larger differences between the nominal and empirical coverages of the STD and DB cubes than models M1 and M2 and rises the changes caused by the inclusion of the parameter uncertainty.In addition, there are no significant changes in empirical coverages for M7 and M8 if the bias correction is implemented, as it happens with M1 and M2 models, for h ¼ 1, and there is evidence of marginal coverages gains when h ¼ 12, being especially visible for highly persistent models.
The main conclusions are similar for T ¼ 250 although, as expected, all empirical coverages get closer to the nominal as Figure 2 shows.We also observe that differences among empirical coverages of bootstrap procedures become less visible for T ¼ 250. 4  According to our Monte Carlo results when forecasting one-step-ahead, the improvement of the bootstrap cubes accuracy is to a large extent due to the parameter uncertainty; see, for example, Clements and Taylor (2001) and Pascual et al. (2004) for a similar result in the context of forecasting univariate AR models.With respect to bias-correction, it seems advantageous when we deal with long horizon forecasts of highly persistent VAR models because the effect of bias raises through powering; see Clements and Taylor (2001), Kim (2004) and Clements and Kim (2007).Moreover, even though bias can be more severe in large dimensional VAR models as shown by Abadir et al.
(1999) and Lawford and Stamatogiannis (2009), its effect on the accuracy of forecast regions seems to be similar to that observed for low-dimensional VAR models.In addition, our Monte Carlo results shed light on the role played by the lag order uncertainty, suggesting that it can be substantial.This, in part, qualifies the arguments by Chatfield (1993); Chatfield (1996) who suggests that model uncertainty could severely worsen the accuracy of forecast regions; see alternatively Kilian (1998a) who found out that the gains of tackling the lag order uncertainty can differ substantially among commonly estimated VAR models, when construing bootstrap forecast intervals for impulse response functions.As expected, given that the LS estimator is consistent, our results confirm that the effects of bias, model, and parameter uncertainty tend to vanish as the sample size increases.
Rows 2-4 of Figures 1 and 2 report the Monte Carlo coverages when the errors are generated by the Student-3, χ 4 , and Bimodal distributions.First, note that, for T ¼ 100, one-step-ahead STD regions do not worsen their coverage accuracy in presence of non-Gaussian errors.Moreover, the one-step-ahead bootstrap cubes that only incorporate the uncertainty about the error distribution do not provide better empirical coverages than STD cubes, regardless of the errors distributions and models, with the exception of the four variate VAR models with Student-3 errors.The last result may be surprising given that the true error distribution is not Gaussian.However, notice that for h ¼ 12, the empirical coverages of the DB cubes are comparable to those of the STD cubes for all errors distributions and models.Therefore, using DB seems like a minimum loss when the errors are other than Gaussian.Also note that, dealing with non-Gaussian errors does not alter the fact that parameter uncertainty is important.Likewise, bias-correction provides only marginal coverages gains when we forecast long horizons.Finally, tackling the lag order uncertainty gives additional accuracy gains for both, short and long horizons, which seem more pronounced when the errors distribution is either χ 4 or Bimodal.Again, Figure 2 shows that, regardless of the errors distributions considered, all empirical coverages improve when we increase the number of observations to T ¼ 250.Nonetheless, DB performs much worse than STD when we consider short horizons and the errors distribution is either χ 4 or Bimodal.
Finally, Figure 3 reports the Monte Carlo results for the 90% Bonferroni cubes when the BIC is used to select the lag order and T ¼ 100.In general, the order of merit of the procedures across models is not affected when the lag order is obtained by BIC (instead of the AIC).For this reason, we do not delve into a detailed description of the 4 When the confidence level is 95%, Monte Carlo results for all the VAR models and procedures are quantitatively the same and, for this reason, are not reported.They are available upon request.
F I G U R E 3 Monte Carlo difference between empirical coverages of 90% h-steps-ahead of bootstrap Bonferroni cubes and nominal coverage, for h ¼ 1 (1st column) and 12 (2nd column).Based on STD (light blue), DB (orange), PB (dark red), BCB (green), and EB (blue) for VAR models with sample sizes T ¼ 100 and Gaussian (1st row), Student-3 (2nd row), χ 4 (3rd row), and Bimodal (4th row) errors.Lag order is selected by using BIC [Colour figure can be viewed at wileyonlinelibrary.com] results. 5Notwithstanding, the comparison of Figures 1  and 3 suggests that the BIC is preferable in terms of empirical coverages for h ¼ 1, whereas AIC is preferred for h ¼ 12.A possible explanation is that the AIC tends to select, on average, a larger lag order than the BIC.Consequently, the AIC uses more in-sample information than the BIC when constructing out-of-sample forecasts.This information may be more valuable when we forecast long instead of short horizons; see, for instance, Marcellino et al. (2006), Clements and Kim (2007), and Müller and Stock (2011).

| EMPIRICAL APPLICATION
A strong belief says that market conditions and economic activity are closely linked.The foundation of this belief is often based on empirical studies that analyze the predictive power of a set of variables representing market conditions on both the monthly returns S&P500 and the growth rate of industrial production; see Fama (1990); Campbell and Diebold (2009); Camilleri et al. (2019), among others.In this section, forecast densities for monthly percent changes of the seasonally adjusted industrial production total index (ipi t Þ, the S&P500 stock market index (sp500 t ), the dividend yield of S&P500 index (div t Þ, and the seasonally adjusted unemployment rate (une t Þ are constructed by using the standard Gaussian approach and the bootstrap procedures described in Section 3. 6 The data is observed from February 1960 to April 2020, accounting for 723 months. 7The parameters are estimated using a fixed rolling windows of T ¼ 100, starting with the estimation period 1960m1-1968m4 and finishing with 2011m1-2019m4.The number of windows used in the estimation is 611.For each window, the lag order is chosen by implementing AIC with an upper bound of 12. Figure 4 plots the lag order and the dominant root for each window.First, it shows that predominant VAR models are of small lag order (the percentage estimated lag order smaller or equal than 2 and 4 is 59.7% and 73.9%, respectively).Second, concerning the dominant root, Figure 5 suggests that the VAR models are highly persistent with a dominant root larger than 0.9 most of the times (98.8%).Thus, in Figure 5, there are many windows for which the dominant root of estimated VAR models is near or equal to one.For instance, the number F I G U R E 4 Estimated lag order (left axis) and dominant root (right axis) of VAR models for ipi t , sp500 t , div t , and une t with fixed windows of T ¼ 100 observations [Colour figure can be viewed at wileyonlinelibrary.com] of estimated VARs with dominant root larger than 0.999 reaches 149 (24.4%).It is worth remarking that bootstrap procedure implemented in this paper is asymptotically valid for a dominant root smaller than one.For this reason, we separately consider the stationary VARs, which we define as those having a dominant root smaller than 0.999, and the non-stationary VAR models, having a dominant root larger than 0.999 or equal to 1. 8 The average dominant root is 0.96 and (roughly) 1 for the stationary and non-stationary VAR models.Finally, we also carried out a multivariate Gaussianity test (Doornik and Hansen (2008)) for the errors for each of the estimation windows and found that Gaussianity is hard to maintained (rejection of Gaussian errors at a significance level of 5% rises to 78.1% of the windows).Furthermore, the presence of non-Gaussian kurtosis is more frequent than non-Gaussian skewness (e.g., multivariate Gaussian kurtosis is rejected at 5% on 65.1% of windows while skewness only on 41.9% of the windows).
Using the h-steps-ahead, for h ¼ 1, …, 12, forecast densities at each window, 90% forecast cubes are obtained.We calculate their empirical coverages by counting the number of observations belonging to them, and we report them in Figure 5 for the stationary and non-stationary VAR models.The upper part of Figure 5 shows that, for stationary VAR models, all procedures' empirical coverages are below the nominal, especially for those dealing only with error uncertainty as STD and DB.As expected, the bootstrap procedures that tackle the parameter uncertainty as PB, BCB, and EB, always get better coverage accuracy than STD and DB, for h ¼ 1 and h ¼ 12.This fact is in line with the results in Figure 1 for model M8, which has four variables, p ¼ 2 and shows high persistence. 9Also, note that PB and BCB achieve roughly the same empirical coverage when forecasting short horizons, although there are marginal gains when forecasting long horizons.Finally, among all procedures that tackle the parameter uncertainty, EB provides the empirical coverages closest to the nominal, regardless of the forecast horizon. 8Beyond the absence of asymptotic validity in the case of non-stationary VAR models, the bootstrap approach can still be useful as an approximation and provide good coverages properties; see, for instance, Kilian (1998a).for the case of impulse response functions. 9Nonetheless, we should be careful because the lag order and the persistence of the estimated VAR model change for each window.
F I G U R E 5 90% empirical coverages of h-steps-ahead Bonferroni cubes, h ¼ 1 left) and 12 (right), for ipi t , sp500 t , div t , and une t with T ¼ 100 observations: STD (light blue), DB (orange), PB (red), BCB (green), and EB (deep blue).Stationary (first row) refers to estimated VAR models with dominant root smaller than 0.999 and non-stationary (second row) to estimated VAR with the largest root greater than 0.999 or equal to 1 [Colour figure can be viewed at wileyonlinelibrary.com]The bottom of Figure 5 reports results for the nonstationary VAR models.Note that the overall effect of uncertainties described for stationary VAR models remains.However, all regions' empirical coverages drop substantially, especially for h ¼ 1, compared with the stationary VAR models.Remarkably, for the non-stationary VARs, Bonferroni regions dealing only with errors uncertainty suffer from a considerable underestimation of coverages.On the other hand, by tackling other sources of uncertainties, we can obtain much better empirical coverages than DB and STD; for instance, see Kim (1999); Kim (2004) for a similar result.
Overall, the results plotted in Figure 5 are in line with those obtained with simulated data.There are large differences when the parameter is included and that the additional gains obtained by tackling bias are not always guaranteed, at least when dealing with small and stationary VAR models.Nonetheless, they become more visible for large and high persistence VAR models when forecasting more than one-step-ahead.Finally, incorporating the lag order uncertainty provides additional empirical coverages gains.

| CONCLUSION
In this paper, we asses the impact of bias, model, and parameter uncertainties on empirical coverages of forecast cubes for several VAR(p) models.With this purpose, we carry out Monte Carlo experiments and construct Bonferroni cubes using bootstrap procedures which are able to differentiate between these uncertainty sources.Our results suggest that a better performance of bootstrap cubes is obtained by considering the parameter and lag order uncertainties.Furthermore, the gains of including the parameter uncertainty are larger when the forecast horizon increases and for large and persistent VAR models, in which the uncertainty surrounding parameters estimates is large; see Riise and Tjøstheim (1984).Likewise, we assess the bias uncertainty and find a lower effect than the one often assumed in applications of bootstrap methods, becoming their gains visible for long-run forecast cubes of persistent VAR models; see, for example, Clements and Kim (2007).Recently, Grabowski et al. (2020) implemented bias-adjustment to higher moments and found additional empirical coverage gains for constructing confidence paths for impulse response coefficients.Nonetheless, it is still an open question how well this approach works when the aim is to forecast a VAR model.Finally, our findings suggest that, when the sample size is large and the parameter uncertainty decreases, it is important to obtain forecast densities that are robust to non-Gaussian errors.The results for highly persistent models are relevant in practice as many macroeconomic time series appear to be closed to a unit root process.

1
Monte Carlo difference between empirical coverages of 90% h-steps-ahead of bootstrap Bonferroni cubes and nominal coverage, for h ¼ 1 (1st column) and 12 (2nd column).Based on STD (light blue), DB (orange), B (dark red), BCB (green), and EB (blue) for VAR models with sample sizes T ¼ 100 and Gaussian (1st row), Student-3 (2nd row), χ 4 (3rd row), and Bimodal (4th row) errors.Lag order is selected by using AIC [Colour figure can be viewed at wileyonlinelibrary.com]F I G U R E 2 Monte Carlo difference between empirical coverages of 90% h-steps-ahead Bonferroni cubes, for h ¼ 1 (1st column) and 12 (2nd column).Based on STD (light blue), DB (orange), PB (dark red/deep mars), BCB (mantle green), and EB (blue) for VAR models with sample sizes T ¼ 250 and Gaussian (1st row), Student-3 (2nd row), χ 4 (3rd row), and Bimodal (4th row) errors.Lag order is selected by using AIC [Colour figure can be viewed at wileyonlinelibrary.com] U R E B . 1 First panel, monthly industrial production total growth rate (ipi t = 100 Â ΔlogðIPI t Þ), second panel, monthly S&P500 returns (sp500 t ¼100 Â ΔlogðS&P500 t ÞÞ, third panel, S&P500 dividend yields growth rate (div t ¼ 100 Â ΔlogðY IELDS t ÞÞ, and fourth panel, unemployment rate (level).Shaded areas correspond to US recession according to NBER A P P END I X B: TIME SERIES PLOT OF SERIES USED IN THE EMPIRICAL APPLICATIONFRESOLI T A B L E 1 VAR(p) models considered in the Monte Carlo experiments