Stationary Bootstrap Based Multi-Step Forecasts for Unrestricted VAR Models

This paper proposes a new asymptotically valid stationary bootstrap procedure to obtain multivariate forecast densities in unrestricted vector autoregressive models. The proposed method is not based on either backward or forward representations, so it can be used for both Gaussian and non-Gaussian models. Also, it is computationally more eﬃcient compared to the available resampling methods. The ﬁnite sample performance of the proposed method is illustrated by extensive Monte Carlo studies as well as a real-data example. Our records reveal that the proposed method is a good competitor or even better than the existing methods based on backward and/or forward representations. optimal lag order is determined by the AIC since it is proposed to determine the best model for forecasting, and the results show that VAR(3) model is the optimal according to AIC.


Introduction
Since the seminal paper of Sims (1980) vector autoregressive (VAR) models have received much attention in econometrics to analyze multivariate time series data. They provide a flexible framework for the analysis of complex dynamics, policy analysis, structural inference and forecasting. Lütkepohl (1991), Watson (1994), and Stock and Watson (2001) are fundamental technical references on this subject and they present an excellent overview of research on VAR models. For the financial applications of VAR model please see Lütkepohl (1991), Campbell et al. (1997), Culbertson (1996), Mills (1999), and Tsay (2001). In this study, we restricted our focus to the forecasting ability of VAR models. VAR model is one of the most commonly used techniques for forecasting, and it has been used in a wide range of applications. See, for example, Batchelor et al. (2007), Gupta et al. (2011), Baumeister and Kilian (2012), D' Agostino et al. (2013), Kilian and Vigfusson (2013), and Hassani et al. (2015) for recent studies on VAR model in econometric applications.
Forecasting with VAR models mostly focuses on the point forecast, which may not provide an accurate inference for a future realization of a time series. As pointed out by Kim (1999), this practice is seriously flawed because the use of point forecasts is meaningless if the extent of the associated uncertainty is unknown. On the other hand, interval forecasts provide better inferences taking into account uncertainty associated with each forecast; see Chatfield (1993), Jore et al. (2010), and Bache et al. (2011). In the context of multivariate time series, the interval forecasts for a given horizon are obtained in the form of ellipsoid and/or Bonferronni cube, see Lütkepohl (1991) and Kim (1999). Traditional approaches to calculate multivariate forecast densities require strong assumptions such as Gaussian forecast errors, known lag order, and model parameters which are generally unknown in practice. Moreover, the constructed forecast densities can be affected due to any departure from the assumptions and may lead us to unreliable results. Such departures from the assumptions are quite possible especially in financial time series. Alternative computational procedures to construct forecast densities without considering distributional assumptions and lag order uncertainties include the use of the well-known resampling methods, such as the bootstrap.
The early studies about obtaining bootstrap forecast densities in VAR models were conducted by Kim (1999) who extends the procedure proposed by Thombs and Schucany (1990) in the context of univariate AR(p) processes. Further, Kim (2001) and Kim (2004) suggested to use bootstrap-after-bootstrap approach to obtain bias-corrected forecast regions. These procedures use backward representation of the model when generating the bootstrap replicates. However, using backward representation increases the complexity of the procedure and restricts the implementation of bootstrap method. To overcome these vexing issues Fresoli et al. (2015) proposed a bootstrap procedure (abbreviated as "FRP") by extending the studies of Pascual et al. (2004Pascual et al. ( , 2006 for constructing multi-step multivariate forecast densities in the context of non-Gaussian unrestricted VAR models. Their procedure avoids the backward representation (it uses forward representation), thus it can be used to obtain multivariate forecast densities in VARMA and VAR-GARCH models.
In this paper, we propose a simple bootstrap procedure, which is originally proposed by Hwang and Shin (2013) for AR(p) models, to obtain multivariate forecast densities in VAR models. It combines the stationary bootstrap of Politis and Romano (1994) and traditional residual based bootstrap method to generate bootstrap replications. The proposed method is not based on either backward or forward representations when generating the bootstrap sample of the original data, hence it requires less computing time compared to such bootstrap methods. Also, it can be applied to VARMA and VAR-GARCH models as in FRP. Our extension works as follows. First, the Ordinary Least Squares (OLS) estimators of the VAR model are obtained and the residuals are calculated. Then, the stationary bootstrap method is applied to the data to obtain the bootstrap sample which is used to calculate the bootstrap estimators of the VAR coefficients. Finally, the future values of the VAR process are obtained by means of bootstrap replicates and quantiles of the Monte Carlo estimates of the generated bootstrap distribution.
The remaining of the paper is organized as follows. In Section 2, we describe the bootstrap procedure based on stationary bootstrap method to construct forecast densities of VAR process. In Section 3, extensive Monte Carlo simulations based on VAR(2) (stationary), VAR(5) (persistent), and VAR(10) (near-cointegrated) models with different parameter values and error distributions are conducted to examine the finite sample performance of the proposed method, and the results are compared with those of FRP. The joint forecast densities of the quarterly US inflation, unemployment, and GDP growth data are obtained, and the results are presented in Section 4. Section 5 concludes the paper.

Methodology
Let us consider an N -dimensional multivariate time series {Y t } = {· · · , Y t−1 , Y t , Y t+1 , · · · }. It follows VAR(p) model with autoregressive order p if where Y t = (Y 1t , · · · , Y N t ) ∈ R N is the available time series consisting of random N × 1 vectors, µ = (µ 1 , · · · , µ N ) is an N × 1 vector of intercept term, Φ i , i = 1, · · · , p, are N × N matrices of coefficients satisfying the stationary condition, and {ε t } is multivariate white noise with nonsingular covariance matrix Σ ε . For j = 1, · · · , N , let There are N linear regression models such that where Φ (j) and ε (j) represent an N p × 1 vector of regression parameters for the jth time series and T vector of errors, respectively. Then, we have a multivariate regression model consisting of N components as follows: The OLS estimator of Φ,Φ = Φ 1 , · · · ,Φ p , is then obtained as follows: Throughout this paper, we obtain the OLS estimates of the VAR(p) parameters taking into account the stationarity restriction of Kilian (1998).
Forecasting forms a large part of multivariate time series analysis. In the context of VAR model, forecasting future observations is similar to forecasting from univariate autoregressive time series models. The point predictor of Y T +h , where h is the forecast horizon, can be obtained by the chain rule of forecasting as follows: If ε t is assumed to be Gaussian then the forecast density of the h-step ahead forecast density may be estimated as follows: The prediction ellipsoid is obtained based on Equation (6). On the other hand, obtaining this ellipsoid becomes complicated when h ≥ 2. Also, the forecast regions obtained by prediction ellipsoid are based on the assumption of Gaussian distributed error terms, so, they may not be reliable when the errors are non-Gaussian. Consequently, a Bonferroni cube was proposed by Lütkepohl (1991) to construct forecast regions. As noted earlier, construction of forecast densities/regions may be affected by a great amount due to any departure from the distributional assumptions and may lead us to unreliable results; see Dufour and Jouini (2006), and Fresoli et al. (2015). Bootstrap-based forecast densities are provide a possible solution to overcome this issue since they are not require the full knowledge of the underlying data and distributional assumptions. In this context, Kim (1999) proposed a bootstrap method based on backward recursion to obtain h-step ahead forecast densities. Kim (2001) and Kim (2004) modified this procedure to obtain bias-corrected bootstrap forecast regions. However, bootstrap methods based on such representation are complicated and time consuming. Also, it limits the applicability of the bootstrap method. Later, Fresoli et al. (2015) proposed a simpler bootstrap method which is based on forward recursion. In this study, we propose a stationary bootstrap based procedure to obtain forecast regions for non-Gaussian unrestricted VAR models. Note that we construct only the Bonferroni cube in our numerical analysis since it provides better forecast regions by taking into account the asymmetry of the distribution compared to forecast ellipsoid.
To start with, consider the observed N -dimensional multivariate time series data , be the blocks of consecutive observations starting from Y i . The observed time series data {Y t } is wrapped around a circle in order to ensure that all starting points have equal probability of selection. Let I 1 , I 2 , · · · be the independetly and identically distributed (iid) discrete uniform random variables on {1, · · · , T } so that P (I 1 = i) = 1/n, for i = 1, · · · , T . Let also L 1 , L 2 , · · · be the iid geometric random variables with parameter ρ such that 0 < ρ < 1 and the probability mass function P (L 1 = ) = ρ(1 − ρ) −1 , for = 1, 2, · · · . We assume that two sets {I 1 , I 2 , · · · } and {L 1 , L 2 , · · · } are independent and ρ → 0 as T ρ → ∞. Then, the stationary bootstrap time series data The stationary bootstrap, among other block based bootstrap methods such as non-overlapping or overlapping block bootstrap, is a useful tool because it has a random block length and captures the stationarity of the original data, see Politis and Romano (1994) for more detailed information. Next, we define our proposed algorithm to construct forecast regions for VAR(p) models.
Step 1. Determine the lag-order p by some information criteria and obtain the OLS estimation of the parameter vector Φ,Φ. Note that we used Akaike information criteria (AIC) to determine the optimal p for all the numerical analysis considered in this paper.
Step 2. Obtain the residualsε t = Y − ZΦ for t = p + 1, · · · , T . LetF ε be the empirical distribution function of the centered and re-scaled residuals.
Step 3. Generate the stationary bootstrap observations Step 4. Calculate the stationary bootstrap estimatorsΦ * = Φ * 1 , · · · ,Φ * p as in (4) but using the bootstrap observations obtained in the previous step, Step 5. Calculate h = 1, 2, · · · step ahead bootstrap replicate of Y T , Y * T +h|T with the following recursion, Step 6. Repeat steps 3-5 B times to obtain bootstrap replicates Ŷ * ,1 T h |T , · · · ,Ŷ * ,B T h |T for each h. Note that B denotes the number of bootstrap replications. Then, the Bonferroni cube with at least (1 − α)100% nominal coverage is obtained as where q * i (γ) is the γth quantile of the generated bootstrap distribution of the ith element of Y * T +h|T .
Remark 1. For a realization of AR(p) process {Y −p+1 , · · · , Y T }, letφ andφ * represent the OLS estimator of the true autoregressive parameter vector φ and its stationary bootstrap counterpart, respectively. If ρ → 0 and T ρ → ∞ as T → ∞ then the Theorem 1 of Hwang and Shin (2013) shows that sup as T → ∞, where p * represents the convergence in conditional probability. Now, consider the realization of N dimensional stationary VAR(p) process, {Y −p+1 , · · · , Y T }, defined as in (1). Then, one can easily show that the stationary bootstrap OLS estimation has the same limiting distribution as the OLS,Φ * p * − →Φ asT → ∞, by extending the above theorem to the VAR(p) case. It can also be shown that, conditional on {Y −p+1 , · · · , Y T }, the h-step ahead stationary bootstrap replicateŶ * T +h|T obtained in Step 5 above converges weakly in probability to Y T +h as T → ∞,Ŷ * T +h|T d − → Y T +h , by using the VAR case of Theorem 3.1 of Thombs and Schucany (1990).

Numerical Results
We conduct a simulation study to investigate the finite sample performance of our proposed stationary bootstrap (SB) forecast regions, and compare our results with the method proposed by Fresoli et al. (2015). This comparison has been made through the coverage probability and volume of the forecast region obtained by the Bonferroni forecast cube. We also calculate the squared errors between the Bonferroni cubes obtained by the bootstrap methods and empirical Bonferroni cube to find out which bootstrap method has the best matched forecast region with the empirical one. It is worth mentioning that we also compared the performance of our proposed method with the bootstrap method based on backward representation proposed by Kim (1999). According to our results, it has similar performance with FRP, and therefore to save space, we only report the comparative study with the proposed method and FRP. The experimental design of our simulation study is similar to those of Fresoli et al. (2015) which is as follows. Three bivariate processes: a stationary VAR(2), a persistent VAR(5), and a near-cointegrated VAR(10); and three error distributions: Gaussian N (0, 1), t(5) (heavy tailed), and χ 2 (4) (skewed); as well as two sample sizes T = 100 and T = 300 (see Appendix for the details) are considered. For each setting M = 2000 Monte Carlo simulations with B = 1999 bootstrap resamples are performed. The significance level α is set to 0.05 to obtain 95% Bonferroni cubes for future h = 1, · · · , 10 observations. The performance metrics are calculated based on the following setting.
• Simulate a bivariate VAR(p) series with the lag order parameters given in Appendix, and generate R = 5000 future values Y T +h for h = 1, · · · , 10. • Calculate bootstrap future valuesŶ * ,b T +h|T for h = 1, · · · , 10 and b = 1, · · · , B. Then, estimate the coverage probabilities (Cp * ) of bootstrap forecast regions forŶ * T +h|T as follows: • Estimate the volume of the bootstrap Bonferroni cube (V Boot ) using the following equation: • Calculate the squared error between the bootstrap and empirical Bonferroni cubes as follows: We note that all the calculations were carried out using R 3.3.2 on an Intel Core i7 6700HQ 2.6GHz PC, and the codes can be obtained from the author upon request.
The results are presented in Figures 1-3 when the sample size T = 100. Our findings show that, regardless of the error distribution, the proposed method have similar coverage performance with FRP when the bivariate data come from the stationary VAR(2) process, and both methods have coverages close to nominal level; see the first row of Figure 1. When the data are generated from the persistent VAR(5) and near-cointegrated VAR(10) processes the coverage performance of the proposed method is better than the FRP especially for long-term forecast horizons as it is shown in the second and third rows of Figure 1. In both cases, although FRP has good coverage values for the Y T +1 its performance gets worse when h ≥ 2. This deterioration may be caused by the fact that the models under consideration are close to the non-stationary bounds especially for the VAR(10) process. The proposed method, on the other hand, is somehow less affected by this distortion. Therefore, our proposed method may be a better choice for long-term estimates compared to FRP. Figures 2 and 3 compare the proposed and FRP methods in terms of the volume of Bonferroni cubes obtained by both methods and their similarities to the empirical Bonferroni cube. Figure 2 shows that, although both methods have almost the same volume when the data is stationary VAR(2) process, the proposed method clearly outperforms FRP when the data follow VAR(5) and/or VAR(10) processes. Also, it is clear from Figure 3 that the Bonferroni cubes obtained by our proposed method are more similar to the empirical Bonferroni cube compared to FRP's Bonferroni cube. In summary, when all the results are interpreted together, we can say that similar and mostly better results can be obtained by the proposed method compared to FRP.
Moreover, we compare the proposed and FRP methods in terms of their computing times, and Figure 4 presents the approximate computing times for various sample sizes based on B = 1999 bootstrap replications and only one Monte Carlo simulation. As it can be seen from Figure 4 that the proposed method has considerably less computing time than FRP. While the computing time of FRP increases with the increasing number of parameters our proposed method has similar computing times for different processes. This is due to fact that FRP uses a recursion to obtain bootstrap resamples while the resamples are obtained directly from the data when using the proposed method. Note that the results for the sample size T = 300 are qualitatively similar to those presented in the paper and therefore not shown to save space. They are available upon request.

Data Example
In this section, we study the forecasting performances of the proposed and FRP methods with a real-world example, quarterly US Real Gross Domestic Product (GDP )-US Inflation rate (ir)-US Civilian Unemployment rate (ur) data which is available at https://www.stlouisfed. org/. The data set, which is consisted a total of 279 observations, observed from 1948Q 1 to 2017Q 3 . The (GDP t ) and (ir t ), respectively, are calculated as GDP t = log(GDP t /GDP t−1 ) and ir t = log(IP D t /IP D t−1 ), where IP D is the Implicit Price Deflator. To obtain out-of-sample forecast cubes for the real data, we divide the full data into the following two parts: the model is constructed based on the observations from 1948Q 1 to 2015Q 1 to calculate 10 steps-ahead forecasts from 2015Q 2 to 2017Q 3 . Table 1  each of the series. The results are as follows: (i) JB test results show that the series are not Gaussian either individually or jointly, (ii) small p-values of ADF test suggest that the series are not stationary processes. Also, we perform the Ljung-Box (LB) test to test for auto-correlations in the original and squared series and smaller p-values indicate that there is a dynamic dependence in the conditional mean. Table 1 also presents the LB statistics of order 10 for the original (Q(10)) and squared series (Q 2 (10)), and they show that the dependence of the second order moments and VAR(10) (third row) models when T = 100 and errors are Gaussian (first column), t(5) (second column) and χ 2 (4) (third column) distributed.
is not as significant as those of generated from the conditional mean dependence since the LB statistics of the squared series are relatively smaller than the ones obtained from the original series. All of our preliminary exploratory analyses suggested to use a VAR model to model the series. The optimal lag order is determined by the AIC since it is proposed to determine the best model for forecasting, and the results show that VAR(3) model is the optimal according to AIC. For each combination of two series, we construct the bootstrap forecast cubes from 2015Q 2 to 2017Q 3 together with the true out-of-sample values which are denoted by a dot, and the results are presented in Figures 5-7. Our findings show that the forecast cubes obtained by both methods are similar and they include all of the observed out-of-sample points. Note that the volumes of the two methods are almost the same for short-term forecast periods while the proposed method produces slightly larger forecast cubes for long-term forecast periods compared to FRP. This is due to fact that the residuals obtained from the proposed method are slightly larger than those of FRP.

Conclusions and Discussion
In this study, we propose an asymptotically valid bootstrap procedure to obtain multivariate forecast densities/regions for future values in vector autoregressive models by extending the stationary bootstrap method proposed by Hwang and Shin (2013). We examine its finite sample performance and compare with the existing method(s) by extensive Monte Carlo studies as well as a real-data example. The prominent feature of the proposed method is that it is not based on either backward or forward representations when generating the bootstrap sample of the original data. Thus, it can be applied for both Gaussian and non-Gaussian models. Also, it requires less computing time compared to existing bootstrap methods based on backward and/or forward representations. Our Monte Carlo experiments show that the similar or even better results can be obtained by the proposed procedure; (i) when the errors are Gaussian and/or the data is clearly stationary the performances of the existing and our method are almost the same; (ii) in case of non Gaussian errors and/or persistent or near cointegrated models, on the other hand, the results produced by the proposed method are better than those of existing models especially for long-term forecast regions, so that our method has better coverage performance with less volume. To summarize, considering its advantages the proposed procedure can be a good alternative to obtain forecast densities in VAR models.