Robust bootstrap prediction intervals for univariate and multivariate autoregressive time series models

Let us consider an $N$ dimensional multivariate time series $\left\{Y_{t}\right\}$. It follows the VAR($p$) model with the AR order $p$ if

where $Y_{t}=\left(Y_{1t},\cdots,Y_{Nt}\right)^{\prime}\in\mathds{R}^{N}$ are the random $N$ vectors of the time series, $\Phi=\left(\Phi_{0},\Phi_{1},\cdots,\Phi_{p}\right)$ is the vector of coefficients, which is a $N\times\left[(p+1)N\right]$ matrix containing the separate $N\times N$ matrices column-wise as blocks, and $\left\{\epsilon_{t}\right\}_{t\geq 1}$ is a normally distributed multivariate white noise with non-singular covariance matrix $\Sigma_{\epsilon}$. Using the lag operator notation, model (2.1) is written as follows:

where $\phi(B)=I-\phi_{1}B-\cdots-\phi_{p}B^{p}$ is matrix a polynomial of order $p$ and $C$ is an $N$ dimensional constant vector.

Following Tsay et al. (2000), a multivariate time series with a multivariate outlier, $\left\{X_{t}\right\}$, can be written as follows:

where $\Delta=\left(\Delta_{1},\cdots,\Delta_{N}\right)^{\top}$ denotes the size of the outlier and $\mathds{1}_{t}^{(s)}$ is the binary indicator variable that characterizes the outlier at time point $s$. The outlier type (AO/IO) is defined based on the matrix polynomial $\eta(B)$. If an AO exists at time $s$, then $\eta(B)=I$. On the other hand, if an IO is present at time $s$, then $\eta(B)=\left(\phi(B)\right)^{-1}$.

An IO that affects innovation $\epsilon_{t}$ can be considered an extraordinary shock influencing the future values of the series $(Y_{t},Y_{t+1},\cdots)$. In contrast, an AO at time point $s$ affects only the corresponding data point $Y_{t}$, not future values of the series, and can have serious effects on the estimated model parameters, residuals, and accuracy of forecasts. Generally, an AO has a much greater effect on forecasting accuracy than an IO, because an IO at time $s$ affects only the corresponding residual $\epsilon_{s}$. Still, an AO, which is far from the bulk of the data, causes a large residual and affects future residuals. For more information about the AOs, IOs, and their effects on the model parameters, consult Fox (1972), Chang et al. (1998) and McQuarrie and Tsai (2003) .

Let us consider the VAR($p$) model given in (2.1). Let $\widehat{Y}_{t}$ and $Y_{T+h}$ for $t=1,\cdots,T$ and $h=1,2,\cdots$ be the fitted values of $Y_{t}$ obtained by a suitable method such as OLS and $h$-step-ahead unobservable future values, respectively. Then, the point forecast of $Y_{T+h}$ is obtained as follows:

where $\widehat{\Phi}=\left(\widehat{\Phi}_{0},\widehat{\Phi}_{1},\cdots,\widehat{\Phi}_{p}\right)$ is the estimated parameter matrix of $\Phi$ and $\widehat{Y}_{T+j|T}=Y_{T+j}$ for $j\leq 0$. Define the forecast error covariance matrix by

where $\Psi_{k}=\sum_{j=1}^{p-1}\Psi_{k-j}\Phi_{j}$ with $\Psi_{0}=\mathbf{I}_{N}$ and $\Phi_{k}=\bm{0}$ for $k>p$ (see, e.g., Lütkepohl, 1991, 2005). For the VAR models, the prediction interval for a given horizon is obtained in the form of an ellipsoid and/or Bonferroni cube (see, e.g., Lütkepohl, 1991; Kim, 1999). Under the assumption of Gaussian distributed errors, the prediction ellipsoid is obtained by the $h$-step-ahead forecast density, which may be estimated as:

where $\widehat{\Sigma}(h)$ is the estimated forecast error covariance matrix. The calculation of this prediction ellipsoid becomes complicated when the forecast horizon $h\geq 2$. Alternatively, a Bonferroni cube, which provides better forecast regions by taking into account the asymmetry of the distribution, is suggested by Lütkepohl (1991). Thus, we calculate only the Bonferroni cube in our numerical analyses.

As noted earlier, this procedure may not be sufficient to obtain accurate prediction intervals since it is based on the normality assumption, conformity to which is generally unknown in practice. Therefore, it may be affected by the departure from normality and lead to unreliable results. Alternatively, the bootstrap method may be used to construct a prediction interval for $Y_{T+h}$ without considering any distributional assumption. The underlying idea behind this method is to approximate the distribution of future observations by drawing resamples conditional on the observed data. Early studies on the construction of bootstrap forecast regions in the VAR models were conducted by Kim (1999), who extended the idea of Thombs and Schucany (1990) to the VAR case. Kim (2001, 2004) proposed a bootstrap-after-bootstrap approach to obtain bias-corrected bootstrap forecast regions. Later, Fresoli et al. (2015) proposed a bootstrap procedure based on the forward recursion to decrease the complexity of the bootstrap method of Kim (1999). In this study, we consider the technique proposed by Fresoli et al. (2015) (abbreviated as “FRP” hereafter), whose algorithm is given below.

• Step 1.

  For a realization of $N$-dimensional multivariate time series $\left\{Y_{t}\right\}_{t=1,\dots,T}$, determine the lag-order parameter $p$ by Akaike Information Criterion (AIC) and calculate the OLS estimator $\widehat{\Phi}$ of the parameter matrix $\Phi$.

• Step 2.

  Compute the residuals $\widehat{\epsilon}_{t}$, for $t=1,\cdots,T$. Let $\widehat{F}_{\epsilon}$ denote the empirical distribution function of the centered and re-scaled residuals, $\widetilde{\widehat{\epsilon}}=\left(\widehat{\epsilon}-\bar{\widehat{\epsilon}}\right)/\hat{\sigma}$, where $\bar{\widehat{\epsilon}}=\frac{1}{T}\sum_{t=1}^{T}\widehat{\epsilon}_{t}$ and $\hat{\sigma}$ is the calculated standard deviation of the residuals, respectively.

• Step 3.

  Generate the bootstrap observations $\left\{Y^{*}_{1},\cdots,Y^{*}_{T}\right\}$ via the following recursion:

  $$Y^{*}_{t}=\widehat{\Phi}_{0}+\sum_{i=1}^{p}\widehat{\Phi}_{i}Y^{*}_{t-i}+\epsilon_{t}^{*},\qquad t=1,\cdots,T,$$

  where $\left\{\epsilon_{t}^{*}\right\}_{t=1}^{T}$ is an independent and identically distributed (iid) sample from $\widehat{F}_{\epsilon}$ and $Y^{*}_{t}=Y_{t}$ for $t=1,\cdots,p$.

• Step 4.

  Compute the OLS estimator $\widehat{\Phi}^{*}=\left(\widehat{\Phi}^{*}_{0},\widehat{\Phi}^{*}_{1},\cdots,\widehat{\Phi}^{*}_{p}\right)$ for the bootstrap observations $\left\{Y^{*}_{1},\cdots,Y^{*}_{T}\right\}$.

• Step 5.

  Conditional on the last $p$ original observations, obtain the $h$-step-ahead bootstrap point forecast of $Y_{T+h}$, given by $Y^{*}_{T+h|T}$ as follows:

  $$Y_{T+h|T}^{*}=\widehat{\Phi}^{*}_{0}+\widehat{\Phi}^{*}_{1}\widehat{Y}^{*}_{T+h-1|T}+\cdots,\widehat{\Phi}^{*}_{p}\widehat{Y}^{*}_{T+h-p|T}+\epsilon^{*}_{T+h},\qquad h\geq 1$$

  where $\widehat{Y}^{*}_{T+h|T}=Y_{T+h}$ for $h\leq 0$ and $\epsilon^{*}_{T+h}$ is randomly drawn from $\widehat{F}_{\epsilon}$.

• Step 6.

  Steps 3-5 are repeated $B$ times to obtain bootstrap replicates $\left\{\widehat{Y}^{*,1}_{T+h|T},\cdots,\widehat{Y}^{*,B}_{T+h|T}\right\}$ for each $h$.

Then, a $100(1-\gamma)\%$ Bonferroni cube is obtained as follows (see, e.g., Kim, 1999)

where $\tau=1-\gamma/(2N)$, $q_{i}^{*}(\tau)$ is the $\tau$^th quantile of the generated bootstrap distribution of the $i$^th element of $\widehat{Y}_{T+h|T}^{*}$, and $\bigtimes$ denotes cartesian product.

For the univariate case, we consider the following AR(1) model:

where $\epsilon_{t}\sim N(0,1)$. The modeling scenarios are adapted such that 1) no clear outlying data points are generated, 2) $C=\left[1\%,5\%,10\%\right]$of only the generated training data are contaminated by deliberately inserted AOs or IOs, and 3) $K=\left[1,2,3\right]$of the generated future values are contaminated by deliberately inserted AOs or IOs. AOs are generated by contaminating $TC$ randomly selected observations with the value 5, while IOs are generated by contaminating $TC$ randomly selected innovations with the value 3. For the first scenario, where no outlying observations are deliberately inserted into the dataset, we consider the parameter values $\phi_{0}=0$ and $\phi_{1}\in\left[-0.9,-0.5,-0.1,0.1,0.5,0.9\right]$ to assess the effect of parameter values on the bootstrap prediction intervals. For the second and third scenarios, the parameter values are fixed at $\phi_{0}=0$ and $\phi_{1}=0.5$, and we check the effects of the number of outliers on the prediction intervals. For these cases, the average coverage probabilities when the bootstrap-based quantiles are replaced with the quantiles of the normal distribution are also calculated to explore whether bootstrapping improves the proposed and other statistical methods. The results are given in Figures 1-6.

Our results show that when no fixed outlying data point is inserted into the dataset: 1) all the bootstrap prediction intervals tend to produce similar average coverage probabilities when the time series data are persistent ($\phi_{1}=\pm 0.9$), with only the Rob-YW slightly differing from other methods; 2) the bootstrap procedure based on the weighted likelihood methodology generally produces improved coverage probabilities than other methods when the generated data are not persistent; 3) when the generated data are not persistent, the proposed method produces coverage probabilities close to the nominal level and generally results in slight under-coverage (over-coverage is observed only for the case when $\phi_{1}=0.5$ and $h>4$, but is observed for all bootstrap prediction intervals based on different methods); and 4) the proposed bootstrap procedure has a smaller prediction interval length than other methods when the data are not persistent, while all bootstrap methods tend to produce similar interval lengths when $\phi_{1}=\pm 0.9$. One exception is that Rob-Reg has smaller interval lengths than other methods when $\phi_{1}=\pm 0.9$ and $h>4$. Note that, when the generated data have no deliberately inserted outliers, the proposed bootstrap method is expected to perform similarly to other methods. However, it performs better than other bootstrap prediction intervals in some cases, because some outlying observations with small magnitudes may be produced by the data generating process. The weighted likelihood method reduces the effects of such observations, giving the proposed method a slightly better performance than other methods in some cases.

When the training data have AOs, over-coverage ($\geq 0.98$, in general) is observed for all bootstrap prediction intervals, excluding the proposed method because of the high forecasting errors produced in line with the outlying points. In contrast, the proposed procedure is less affected by these points and provides coverage probabilities close to the nominal level. In addition, the proposed idea has significantly narrower prediction intervals than other methods when AOs are present in the data (see Figure 3).

When the generated training data have IOs, compared with other methods, the proposed bootstrap method generally performs poorly for small $h$ in terms of empirical coverage probability. Under-coverage is observed for the proposed method while the other methods result in over-coverage. As in the AO case, the proposed bootstrap method produces significantly narrower prediction interval lengths compared with other methods. When comparing the performance of the bootstrap methods for the AO and IO cases (Figures 3 and 4), the results demonstrate that the bootstrap methods produce slightly lower coverage values when AOs are present in the generated datasets than those of IOs. This is due to fact that the bootstrap methods tend to produce wider prediction intervals when the data are contaminated by AOs compared to IO case. As is mentioned in Section 2.1, an AO at time $t$ influences corresponding residual as well as the future residuals (at time $t+1,\cdots$) while an IO at time $t$ influences only the corresponding residual. Therefore, compared with the IO case, more high-valued forecast errors are present when the data have AOs. As a result, the bootstrap methods produces wider prediction intervals when AOs are present in the data compared with IOs, and thus, the prediction intervals cover more future values.

Moreover, our findings for the coverage probabilities indicate that bootstrap-based quantiles generally perform better than the quantiles of the normal distribution when the time series is generated without outliers. When the training data are contaminated by AOs, the proposed bootstrap method outperforms normal quantiles for small forecast horizons, while the quantiles obtained from the proposed method and the normal distribution produce similar coverage performances for mid-term and long-term forecast horizons. In contrast, when the training data are contaminated by IOs, the quantiles of the normal distribution perform better than all bootstrap-based quantiles for moderate and large values of $h$, while the bootstrap-based quantiles perform better when $h$ is small. Our conclusions about the comparison between bootstrap-based and normal quantiles are restricted to data structures considered in this study, and may not be readily generalizable.

When the generated future values include AOs (see Figure 5), under-coverage is observed for all the bootstrap methods and normal quantiles. In this case, all the methods perform considerably less in terms of coverage probability compared with the second scenario, and their coverage performances get worse as the number of outliers increase. On the other hand, when the generated future observations have IOs (see Figure 6), for all bootstrap methods, under-coverage is observed for short-term forecast horizons while over-coverage is observed for long-term forecast horizons. For this scenario, the proposed bootstrap method produces smaller prediction interval length, among others.

The experimental design of our simulation study for the VAR case is somewhat similar to that of Fresoli et al. (2015). We consider the following VAR(2) model:

where

The distribution of the error process depends on the modeling scenarios; three error distributions-Gaussian $N(0,1)$, $t_{5}$ (heavy-tailed) and $\chi^{2}(4)$ (skewed)-are considered when the data have no clear outliers while the errors are distributed as $N(0,1)$ and $t_{5}$ when AOs or IOs are present in the data. To generate AOs, $C=\left[1\%,5\%,10\%\right]$ of the generated data are randomly selected, and both components contaminated by the value 5 while $TC$ innovation terms are randomly selected and only the first component of the innovations are contaminated by the value 5 to generate the IOs.

Our results are reported in Figures 7-11. When no outlier is present in the data, all bootstrap methods tend to produce similar coverage probabilities and volume of forecast cubes, as shown by Figure 7. Only the Rob-VAR produces slightly different coverage probabilities compared with the proposed and OLS methods when the errors follow $N(0,1)$ and $t_{5}$ distributions. Figure 7 also shows that the proposed method has fewer errors than the other two methods when the errors are generated from $N(0,1)$ and $t_{5}$ distributions. However, it produces larger errors than the other methods when the errors follow a $\chi^{2}(4)$ distribution.

When the data have AOs/IOs and the error terms follow $N(0,1)$ distribution, the proposed method has significantly better performance compared with the other two methods. In other words, the proposed method produces better coverage probabilities, less volume, and Bonferroni forecast cubes more similar to the empirical Bonferroni forecast cube compared with other bootstrap methods. All the methods result in under-coverage when outliers contaminate the data; however, the coverage probabilities produced by the proposed method are much closer to the nominal level than those of other bootstrap methods (see Figures 8 and 9). When the data have AOs/IOs but $t_{5}$ distributed error terms (Figures 10 and 11), the proposed method still has a superior performance over other methods, but it performs slightly worse compared with the case where the errors follow $N(0,1)$ distribution.

Original monthly Brazil/U.S. dollar exchange rate data were obtained starting from February 1, 1999 to August 1, 2018 (235 observations), and a logarithmic return series generated via

where $y_{t}$ represents the observed exchange rate at the $t$^th month. Time series plots of the original and logarithmically transformed exchange rate data are presented in Appendix, which shows several clear outlying observations.

The $p$-values obtained by the Ljung-Box (LB) and augmented Dickey-Fuller (ADF) $t$-statistics tests ($p$-value $\leq 0.01$) suggest that the return series is a stationary process with mean zero. The autocorrelation and partial-autocorrelation plots of the return series (not presented here) show that an AR model may be appropriate to model the return series. The optimal lag order $p$ for this series is obtained as 3, AR(3), via the AIC. In addition, the structural break test based on $F$ statistics in the R package strucchange of Zeileis et al. (2002) is applied to the returns series to check if there is a structural break in the series. In this test, a $F$ statistic is computed for every predetermined potential change point and the structural break test is performed using the calculated $F$ statistics. For the return series of monthly Brazil/U.S. dollar exchange rate data, this test procedure produce the $p$-value$=0.188$, which indicates that there is no clear structural break in the data.

To obtain out-of-sample bootstrap prediction intervals for the return series, we divide the data into two parts; we estimate the model based on observations from February 1, 1999, to August 1, 2017, to construct $12$-step-ahead prediction intervals from September 1, 2017, to August 1, 2018. Our results are shown in Figure 12, which shows the prediction interval produced by the proposed method is narrower than those obtained by other bootstrap methods, and all prediction intervals include all the observed future values.

Quarterly GDP and IR data for Turkey, which consist of 76 observations, were obtained for 1999$Q_{1}$ to 2017$Q_{3}$. The GDP and IR data were then logarithmically transformed as $GDP_{t}=\ln(GDP_{t}/GDP_{t-1})$ and $IR_{t}=\ln(IR_{t}/IR_{t-1})$, respectively. The time series plots of the original and transformed data are given in Appendix, which shows several outliers for both series.

For this dataset, the model is constructed based on data spanning from 1999$Q_{1}$ to 2016$Q_{3}$ to calculate Bonferroni cubes for four-step-ahead forecasts from 2016$Q_{4}$ to 2017$Q_{3}$. Table 1 presents the sample statistics and Jarque-Bera (JB) normality and ADF stationarity test results. The results indicate that: 1) the data do not follow a Gaussian process, either individually or jointly, and 2) both GDP and IR series are stationary processes according to ADF test results. Moreover, the results of the LB test (not presented here) indicate that there is a dynamic dependence on the conditional mean for this dataset. The optimal lag-order is determined by AIC, with the results showing that a VAR(3) model is optimal.

The constructed Bonferroni forecast cubes and observed out-of-sample values are presented in Figure 13. The bootstrap forecast cubes have considerably less volume than the Bonferroni cubes obtained by other bootstrap methods that use the OLS and multivariate least trimmed squares. Because of the high forecast errors caused by the presence of outliers, the weighted likelihood-based bootstrap method is less affected by the outliers and it thus provides Bonferroni cubes with less volume than the other two bootstrap methods.