Time-varying NoVaS Versus GARCH: Point Prediction, Volatility Estimation and Prediction Intervals

2.1.1 The Benchmark: GARCH(1,1)

Let Y₁, …, Y_n be an observed stretch from a financial returns time series {Y_t, t ∈ Z}. To start with, assume that {Y_t} is (strictly) stationary with mean zero, which implies that trends and other non-stationarities have been successfully removed.

The standard ARCH(p) model is described by an equation of the type:

where the series {Z_t} is assumed to be independent and identically distributed (i.i.d.), and p is a positive integer indicating the order of the model. Engle (1982) assumed that the errors Z_t have a N(0, 1) distribution; later on other researchers considered instead a t–distribution (standardized to unit variance) for the errors—see e.g. Shephard (1996).

Let ℱn be short-hand for the observed information set, i.e., ℱn={Yt, 1≤t≤n}. Under the above ARCH(p) model, the L₂ optimal predictor of Yn+12 based on ℱn is

This conditional expectation E(Yn+12|ℱn) is commonly referred to as the volatility. although the same term is sometimes also used for its square root.

A standard GARCH(1,1) model is described by the equation:

where the series {Z_t} is i.i.d. (0,1), and the parameters A, B, C are assumed nonnegative. The quantity ht2=E(Yt2|ℱt−1) is the volatility of the GARCH model; it is also the L₂ optimal predictor of Yt2 given ℱt−1.

Back-solving in the right-hand-side of Eq. (3), it is easy to show that the GARCH(1,1) model is tantamount to the ARCH model of Eq. (1) with p = ∞ and the following identifications:

Under the objective of L₁ optimal prediction, the optimal predictor is the conditional median—not the conditional expectation. For an ARCH(p) process, the L₁ optimal predictor of Yn+12 given ℱn is

Under the usual causality assumption, the error Z_t is independent of ℱt−1 for all t. Hence, Median(Zn+12|ℱn)=Median(Zn+12).

Furthermore, the aforementioned equivalence of GARCH(1,1) to a particular ARCH(∞) model implies that Eq. (5) would also give the L₁ optimal GARCH(1,1) predictor of Yn+12 by allowing p = ∞ with the ARCH coefficients a, a₁, a₂, … having the structure given by Eq. (4). Alternatively, under the GARCH(1,1) model we can simply write

Note that if Z_t has an (approximately) symmetric unimodal distribution, then Zt2 is a positively skewed random variable. It follows that Median(Zn+12)≤E(Zn+12)=1. For example, Median(Zn+12)≃0.45 if Z_t is N(0, 1), whereas Median(Zn+12)≃0.53 if Z_t has a t distribution with 5 degrees of freedom. Hence, in both ARCH and GARCH cases, the L₁ optimal predictor of Yn+12 is quite smaller than its L₂ analog.

The GARCH(1,1) is by far the most popular of the ARCH/GARCH models, and typically forms the benchmark for modeling financial returns. That is why we compare the predictive ability of NoVaS methodology to that of GARCH(1,1) in what follows.

Remark. If the objective is point prediction of Yt2 given ℱt−1, should one choose and L₂ and L₁ optimal predictor? The crux here is that financial returns are governed by a fat-tailed distribution; it is generally accepted that they have a finite variance but it is not clear that they a finite fourth moment. An L₂ measure of prediction of the squares is essentially a fourth moment; if it is infinite, it does not make sense to try to minimize it! The paper by Politis (2007) brought attention to this problem and adopting an L₁ measure of prediction of the squares, i.e. prediction by conditional median (not mean); see also Ch. 10.4.1 of Politis (2015).

2.1.2 NoVaS Methodology

Let us continue considering a zero mean and (strictly) stationary financial return time series {Y_t, t ∈ Z}. The NoVaS methodology attempts to map the dataset Y₁, …, Y_n to a sample from a Gaussian stationary series.

The starting point is the ARCH model of Eq. (1) that can be re-written as

Engle (1982) assumed that the errors Z_t are i.i.d.N(0, 1); this assumption was later challenged by empirical researchers but we re-visit it here. Note that the ratio given in Eq. (7) can be interpreted as an attempt to “studentize” the return Y_t by dividing with a time-localized measure of the standard deviation of Y_t. However, there seems to be no reason to exclude the value of Y_t from an empirical, causal estimate of the standard deviation of Y_t; recall that a causal estimate is one involving present and past data only, i.e., the data {Y_s, s ≤ t}.

Under this rationale, Politis (2003) defined a new “studentized” quantity as follows:

In the above, st−12 is an estimator of the unconditional variance σY2=Var(Yt) based on the data up to (but not including^[1]) time t; under the zero mean assumption for Y_t, the natural estimator is st−12=(t−1)−1∑k=1t−1Yk2.

The definition in Eq. (8) describes the basic normalizing and variance-stabilizing transformation (NoVaS) under which the data series {Y_t} is mapped to the new series {W_t,a}. The order p(≥0) and the vector of nonnegative parameters (α, a₀, …, a_p) are chosen by the practitioner with the twin goals of normalization and variance stabilization.

The NoVaS transformation of Eq. (8) can be re-arranged to yield:

The only essential difference between the NoVaS of Eq. (9) and the ARCH of Eq. (1) is the presence of the term Yt2 in the right-hand-side paired with the coefficient a₀. Replacing the term a in Eq. (1) by the term αst−12 in Eq. (9) is natural since the former has—by necessity—units of variance; in other words, the term a in Eq. (1) is not scale invariant, whereas the term α in Eq. (9) is.

The target of variance stabilization for the new series {W_t,a} relies on the construction of a local estimator of scale for studentization, and requires the “unbiasedness” condition:

Eq. (10) has the interesting implication that the {W_t,a} series can be assumed to have an (unconditional) variance that is (approximately) unity. Nevertheless, note that p and α,  a₀,  …,  a_p must be carefully chosen to achieve a degree of conditional homoscedasticity as well; to do this, one must necessarily take p small enough—as well as α small enough or even equal to zero—so that a local (as opposed to global) estimator of scale is obtained.

Politis (2003) provided two structures for the a_i coefficients satisfying Eq. (10). One is to let α = 0 and ai=1/(p+1) for all 0 ≤ i ≤ p; this specification is called the simple NoVaS transformation. The other one is given by the exponential decay NoVaS where α = 0 and ai=c′e−ci for all 0 ≤ i ≤ p. The simple and exponential NoVaS schemes are very intuitive as they correspond to the two popular time series methods of obtaining a “local”, one-sided average, namely a moving average (of the last p+1 values) and “exponential smoothing”; see e.g. Hamilton (1994).

With regards the target of normalizing the distribution of the series {W_t,a}, note that

since all the parameters are nonnegative; therefore, the random variable W_t,a is bounded:

So one must be careful to ensure that the {W_t,a} variables have a large enough range such that the boundedness is not seen as spoiling the normality. Thus, we also require

for some appropriate C of the practitioner’s choice. Recalling that 99.7% of the mass of the N(0,1) distribution is found in the range ±3, the simple choice C = 3 can be suggested; this choice seems to work reasonably well for the usual samples sizes.

2.1.3 Basic NoVaS Schemes

We then proceed to choose p and α, a₀, a₁, …, a_p (the parameters needed to identify) with the optimization goal of making the {W_t,a} series as close to normal as possible. To quantify this goal, one can attempt to minimize a (pseudo)distance measuring departure of the transformed data from normality. Recall that it is a matter of common practice to assume that the distribution of financial returns is symmetric (at least to a first approximation) and, therefore, the skewness of financial returns is often ignored. In contrast, the kurtosis is typically quite large, indicating a heavy tailed distribution. Hence, the kurtosis can serve as a simple measure of the departure of a (non-skewed) dataset from normality. Other measures of non-normality can be used, e.g., the Shapiro–Wilk statistic or Kolmogorov distance; see Politis and Thomakos (2008, 2013). In practice, all these measures give similar results; hence, we focus on the kurtosis measure which is also easy to interpret.

Let KURT_n(Y) denote the empirical kurtosis of a dataset {Y_t, t = 1, …, n}, i.e.,

where Y¯=n−1∑t=1nYt denotes the sample mean.

The following two algorithms can be used to select the optimal parameters for NoVaS; see Politis (2003). Note that the only free parameter in Simple NoVaS is the order p; therefore, the Simple NoVaS transformation will be denoted by Wt,pS. In the Exponential NoVaS, to specify all the a_is, one just needs to specify the two parameters p and c > 0, in view of Eq. (10). However, because of the exponential decay, the parameter p is now of secondary significance; thus, we concisely denote the exponential NoVaS transformation by Wt,cE.

2.1.4 Generalized NoVaS schemes

Simple NoVaS has only one parameter, and the same is (effectively) true for Exponential NoVaS. Although many different multi-parameter NoVaS schemes can be devised, we now elaborate on the possibility of a nonzero value for the parameter α (say α₁, α₂, …, α_K, that span a subset of the interval [0,1]) in Eq. (8) in connection with the Simple and Exponential NoVaS. We thus define the Generalized Simple (GS) and Generalized Exponential (GE) NoVaS denoted by Wt;p,αGS and Wt;c,αGE indicating their respective two free parameters; both are based on Eq. (8).

2.1.5 NoVaS Prediction Equation

Suppose that the NoVaS parameters, i.e., the order p(≥0) and the parameters α,  a₀,…,  a_p have already been chosen. Re-arrange the NoVaS Eq. (8) to yield:

and

Let g(⋅) be some (measurable) function of interest; examples include g₀(x) = x, g1(x)=|x|, and g2(x)=x2, the latter being the function of interest for volatility prediction. From Eq.(14) it follows that the predictive (given ℱn) distribution of g(Yn+1) is identical to the distribution of the random variable

where An=αsn2+∑i=1paiYn+1−i2 is treated as a constant given the past ℱn, and the random variable W has the same distribution as the conditional (on ℱn) distribution of the random variable Wn+1,a.

Therefore, the L₁ optimal prediction of g (Y_n+1) given ℱn is given by the median of the conditional (given ℱn) distribution of g(Y_n+1), i.e.,

Specializing to the case of our interest, i.e., volatility prediction and the function g2(x)=x2 yields the NoVaS predictor:

Remark. By construction, the series {W_t,a} is approximately Gaussian but can be correlated. Nevertheless, applied work involving numerous real datasets has invariably shown that the series {W_t,a} is effectively white noise, and hence it is i.i.d. due to Gaussianity; see Politis (2015) and the references therein. As a result, the conditioning on ℱn can be dropped in both Eqs. (16) and (17). In this case, a quantity such as μ₂ of Eq. (17) is the unconditional median, and can readily by estimated by the sample median of {W_t,a for t = p+1, …, n}.

The above remark helps clarify the philosophical distinction between model-based and model-free inference. For example, defining Ut,a=Wt,a/1−a0Wt,a2 we can re-write Eq. (14) as

which formally looks like the defining equation of an ARCH(p) model. The difference is that in the model-based approach, one assumes that the U_t,a are i.i.d. in an equation such as (18), and then proceeds to estimate the coefficients a_i by (say) Maximum Likelihood. By contrast, in the model-free setting of NoVaS, the W_t,a (and therefore the U_t,a as well) are i.i.d. by construction, not by assumption. Furthermore, this construction has already taken place before being able to even write down Eq. (18). In other words, the predictive Eq. (18) is the conclusion of the model-free fitting of NoVaS, where in a model-based setup such an equation would be the starting point.

2.2.1 Time-varying GARCH and Time-varying NoVaS

Intuitively, a locally stationary process can be considered as approximately stationary over a short time window which has led to the notion of Time-varying ARCH/GARCH models. In particular, the theory of Time-varying ARCH (TV-ARCH) processes was developed, and the asymptotic properties of weighted quasi-likelihood estimators were studied in detail by Dahlhaus and Rao (2006).

The analysis of a Time-varying ARCH/GARCH model can be based on the premise of local stationarity. For example, in order to predict g(Y_t+1) based on ℱt via a Time–varying GARCH(1,1) model, we can simply fit the GARCH(1,1) model of Eq. (3) using the “windowed” data Yt−b+1, …, Yt; as a result, the coefficients of the fitted GARCH(1,1) model are varying with time. Here, the window size b should be large enough so that accurate estimation of the GARCH parameters is possible based on the subseries Yt−b+1, …, Yt but small enough so that such a subseries can plausibly be considered stationary.

In a similar vein, we can predict g(Y_t+1) by fitting one of the NoVaS algorithms (Simple versus Exponential, Generalized or not) just using the “windowed” data Yt−b+1, …, Yt. In so doing, we are constructing a Time-varying NoVaS (TV-NoVaS) transformation as discussed in Politis (2015). Recall that in extensive numerical work, Politis and Thomakos (2008, 2013) showed that NoVaS fitting can be done more efficiently than GARCH fitting by (numerical) MLE. Thus, it is conjectured that TV-NoVaS may be able to capture a slowly changing stochastic structure in a more flexible manner; stated in different term, the window size b required for accurate NoVaS fitting could be smaller than the one required for accurate GARCH fitting, thus allowing us to capture the time-varying stochastic structure more accurately.

2.2.2 Structural Breaks and Change Points

A different form of non-stationarity is due to the possible presence of structural breaks, i.e., change points, occurring at some isolated time points. Kokoszka and Leipus (2000) and Berkes, Horváth, and Kokoszka (2004) have studied the detection/estimation of change points in ARCH/GARCH modeling. Mikosch and Starica (2004) and Stărică and Granger (2005) show the interesting effects that an undetected change point may have on our interpretation and analysis of ARCH/GARCH modeling. Polzehl and Spokoiny (2006) show that time-varying models can give a better predictive power for time series with change points.

Fitting ARCH/GARCH models with change points can be quite cumbersome from a practical point of view. So it is of interest to investigate the performance of the aforementioned time-varying GARCH and NoVaS methods given data that are stationary except for breaks. If data over short windows are used, then TV-GARCH and TV-NoVaS should perform generally well except in the immediate neighborhood of the change points. Hence, in the simulations that follow, we also include a structural break model in order to see the effect of an undetected change point on the performance of TV-GARCH and TV-NoVaS predictors for squared returns.

2.3.1 Simulation Design

For the simulation, 500 datasets Y̲n=(Y1,…,Yn)′ were constructed using either a TV-GARCH or a change point GARCH (CP-GARCH) model; these were defined using the standard GARCH model of Eq. (3) as building block with C = 10⁻⁵. The i.i.d. errors Z_t associated with a GARCH model are often assumed to have a Student t₅ distribution. Instead, we use the simple assignment Zt∼ i.i.d. N(0, 1) in the simulation in order to give the GARCH fitting a “head-start”; using a GARCH with normal errors as DGP should facilitate the convergence of the numerical (Gaussian) MLE employed in fitting the TV-GARCH model.

The exact models used to generate simulated data are as follows:

1. CP-GARCH: For t≤n/2, let A = 0.10 and B = 0.73; for t>n/2, let A = 0.05 and B = 0.93. These values are close to the ones used by Mikosch and Starica (2004).

2. TV-GARCH: The value of A decreases as a linear function of t, starting at A = 0.10 for t = 1, and ending at A = 0.05 for t = n. At the same time, the value of B increases as a linear function of t, starting at B = 0.73 for t = 1, and ending at B = 0.93 for t = n.

The difference between the CP-GARCH and the TV-GARCH model is an abrupt versus a smooth transition of the parameters spanning the same values. Some more information on the simulation follows.

1. The prediction method employed was the conditional median obtained either from a TV-GARCH model with normal errors fitted by windowed Gaussian MLE, or via TV-NoVaS (Simple or Exponential); in either case, two window sizes were tried out, namely b = 125 or 250.

2. The sample size was n = 1001 corresponding to about 4 years of daily data; so the choices b = 125 and 250 correspond to 6 and 12 months respectively.

3. Training period for all methods was 250, i.e., the experiment amounted to predicting Yt+12 from the “windowed” data Yt−b+1, …, Yt for t = 250,  251,  … ,  1000.

4. Updating (re-estimation) of all methods would ideally be for each t = 250,  251,  …,  1000. To save computing time, updating in the simulation was only performed for t being an integer multiple of 50. In fairness, the performance of predictions was recorded and compared only at the moment of updating the model, i.e., at time points 250,  300,  350, … ,  1000. For each point, we give the mean absolute deviation (MAD) of the prediction error at the update time point averaged over the 500 replications.

2.3.2 Results and Conclusions

Figure 1 shows the MAD of volatility prediction errors of TV-GARCH as compared to that of TV-NoVaS (Simple or Exponential) with data from model CP-GARCH for the 16 time points where the updating and prediction occurred, i.e., the time points 250,  300,  350, …,  1000. The left panel depicts the case b = 125 while the right panel depicts the case b = 250. Figure 2 is similar but using data generated by a TV-GARCH model instead.

Figure 1:

MAD of prediction of squared returns obtained by fitting TV-GARCH versus TV-NoVaS; data from CP-GARCH model.

Figure 2:

MAD of prediction of squared returns obtained by fitting TV-GARCH versus TV-NoVaS; data from TV-GARCH model.

Some conclusions are as follows:

1. The CP-GARCH model of Figure 1 is stationary up to time point 500 where the break occurs. Hence, time points 250, 300, 350, 400, 450, and 500 in the left panel of Figure 1 corroborate the aforementioned fact that NoVaS (Simple or Exponential) beats GARCH for prediction of squared returns even if the data generating model is (stationary) GARCH as long as the sample size available for model-fitting is small—equal to 125 in this case which is the size of the local window. The corresponding points in the right panel of Figure 1 indicate that GARCH manages to do as well as (or better than)^[3] NoVaS when the effective sample size is increased to 250.

2. Figure 1 shows that the change point at t = 500 wreaks havoc in TV-GARCH model fitting and the associated predictions; this adds another dimension to the observations of Mikosch and Starica (2004). By contrast, both NoVaS methods seem to adapt immediately to the new regime that occurs after the unknown/undetected change point.

3. Figure 2 shows that TV-NoVaS (Simple or Exponential) beats TV-GARCH for prediction of squared returns even when the data generating model is indeed TV-GARCH. Not only is the MAD of prediction of TV-NoVaS just a small fraction of that of TV-GARCH, but the wild swings associated with the latter indicate the inherent instability of GARCH model-fitting; this instability is prominent even in this simplistic case where the errors have a true Gaussian distribution, and Gaussian MLE is used for estimating just the three GARCH parameters.

4. As seen in both Figures 1 and 2, the performance of Simple NoVaS is practically indistinguishable from that of Exponential NoVaS although upon closer look the latter appears to be marginally better.

5. It is possible to also consider Time-varying version of Generalized NoVaS both Simple and Exponential. However, in the context of the present simulation, the results of TV GS NoVaS and TV GE NoVaS are practically indistinguishable from the results of the Simple and Exp-Novas; they were not added to Figures 1 and 2 in order not to clutter the pictures and maintain their readability.

Remark. GARCH model fitting via (numerical) MLE is well-known to be problematic unless the sample size is quite large, e.g. more than 300. However, as discussed in the Remark given at the end of Section 2.1.3, NoVaS fitting is straightforward and does not employ numerical MLE. This is what makes NoVaS outperform GARCH even when the model is GARCH with a moderate sample size; see the left panel of Figure 1 where for the first 500 entries the model is pure GARCH. In this sense, is it expected that TV-GARCH is not performing well with locally stationary data (see Figure 2), since the MLE fitting has to be done over short windows.

4.3.1 Three Illustrative Datasets

In the work of prediction intervals, we also focus on another three real world datasets of daily returns taken from a foreign exchange rate, a stock price, and a stock index. A brief description of these datasets is as follows; for more details, see Politis (2015).

1. Example 1: Foreign exchange rate. Daily returns from the Yen versus Dollar exchange rate from January 1, 1988 to August 1, 2002; the data were downloaded from Datastream. The sample size is n = 3600 (weekends and holidays are excluded).

2. Example 2: Stock index. Daily returns of the S&P500 stock index from October 1, 1983 to August 30, 1991; the data are available as part of the GARCH module in Splus. The sample size is n = 2000.

3. Example 3: Stock price. Daily returns of the IBM stock price from February 1, 1984 to December 31, 1991; the data are again available as part of the GARCH module in Splus. The sample size is n = 2000.

4.3.2 Simulation

In the following simulation work, we will compare the performance in interval prediction of squared returns by using Simple-NoVaS, Exp-NoVaS, Limit Model-Free with Simple-NoVaS method (LMF Simple-NoVaS), Limit Model-Free Exp-NoVaS method (LMF Exp-NoVaS), Generalized Simple-NoVaS method (GS-NoVaS), Generalized Exp-NoVaS method (GE-NoVaS), Limit Model-Free with Generalized Simple-NoVaS method (LMF GS-NoVaS) and Limit Model-Free with Generalized Exp-NoVaS method (LMF GE-NoVaS) method with that by using GARCH(1,1) models with normal errors. Each interval is constructed based on the windowed data series. So all the methods and models here are still time-varying models. For the window size, we will use 125 and 250, as are worked in the point predictions. For prediction intervals using a GARCH(1,1) model, we use the Algorithm 9.2.1 of Model-Based prediction intervals for Y_n+1 from Politis (2015); see also Pan and Politis (2016).

As in Section 2, we employ the same TV-GARCH(1,1) and CP-GARCH(1,1) models to generate the simulated datasets. For consistency, two stationary processes are generated by the following two standard GARCH(1,1) models:

• Model 1. Yt=σtϵt, σt2=.00001+.93σt−12+.05Yt−12, {ϵt}∼i.i.d. N(0, 1).

• Model 2. Yt=σtϵt, σt2=.00001+.73σt−12+.10Yt−12, {ϵt}∼i.i.d. N(0, 1).

Each simulated dataset has size n = 1000, and the window sizes are b = 125 and 250, similar to those used in Section 2 of point predictions. For the three real world data discussed in Section 4.3.1, we use window sizes b = 250 and 500 because their full sample size is bigger. For computational reasons, we chose B = 500 for bootstrap resampling throughout; a practitioner that has to deal with a single dataset could well afford to use a larger B, maybe even 10 times larger.

For each dataset, n−b windowed datasets of size b are extracted. For each windowed dataset, one of the bootstrap methods—based on GARCH(1,1) or on different NoVaS algorithms—was used to create B bootstrap sample paths and B one-step ahead “future” values denoted by Y(b+1,j) for j = 1,  2,  … ,  B; The bootstrap prediction interval (L_i, U_i) was constructed for the “future” value Y_(b+1) of windowed dataset i, here i = b+1, b+2, …, n. The corresponding empirical average coverage level (CVR) and the average length (LEN) of the constructed intervals, and the standard error (St.err) associated with each length of the constructed intervals are calculated as

where

4.3.3 Discussion of Results

Our results are summarized in Tables 3 and 4 for simulated stationary datasets, in Tables 5 and 6 for simulated datasets generated by CP-GARCH(1,1) and TV-GARCH(1,1) models respectively and in Tables 7–9 for the three real world datasets. Each table has two subtables with different window sizes. The first two lines of each subtable are the results using the Simple NoVaS (Simple-NoVaS) and Exponential NoVaS (Exp-NoVaS) methods. Lines 3 and 4 of each subtable are the results using Limit Model-Free with Simple-NoVaS (LMF Simple-NoVaS) and Limit Model-Free Exp-NoVaS (LMF Exp-NoVaS) methods. Lines 5 and 6 of each subtable are the results of Generalized Simple-NoVaS (GS-NoVaS) and Generalized Exp-NoVaS (GE-NoVaS) methods. Lines 7 and 8 of each subtable are for Limit Model-Free with Generalized Simple-NoVaS method (LMF GS-NoVaS) and Limit Model-Free with Generalized Exp-NoVaS (LMF GE-NoVaS) method. The last line of each subtable gives the results by using GARCH(1,1) models with normal errors.

For simulated stationary time series data by GARCH(1,1) models, Tables 3 and 4 show that the performance of different NoVaS methods is better than that of GARCH(1,1), that is, NoVaS methods have average coverages very close to or equal to the nominal ones. Remarkably, the Generalized Simple and/or Exponential NoVaS methods beats Simple/Exp-NoVaS as well as LMF Simple/Exp-NoVaS. In addition, there are not significant differences between Limit Model-Free NoVaS and Generalized NoVaS methods when the data is stationary. It is also shown that the window size has little effect—the coverages are still close to nominal ones—for the performance of NoVaS methodology, while the reaction of GARCH(1,1) to the change of the window size is relatively large. We can see that the average coverages get much closer to the nominal ones as the window size is increased from 125 to 250.

Tables 6 and 7 give the results for data generated from time-varying GARCH(1,1) and CP-GARCH(1,1) models with normal errors respectively. It is shown that in each table, all NoVaS methods have the average coverages much closer to the nominal ones than GARCH(1,1) in both subtables for each table. When the window size is increased to 250, GARCH(1,1) models perform better as expected from earlier work of Politis and Thomakos (2008, 2013). Meanwhile, the performance of NoVaS methods is not sensitive to changes in the sample size; this finding is in accordance with the results for point prediction of squared returns of Section 2. However, this effect of window size on GARCH(1,1) under data with nonstationarity is bigger than when the data are stationarity; compare to Tables 3 and 4.

It is also worthwhile to note that Limit Model-Free NoVaS and Generalized NoVaS methods can sometimes capture the nominal coverages exactly with smaller average lengths and standard deviations for the predicted intervals. Also, we can find that Limit Model-Free NoVaS methods perform better than Simple-NoVaS and Exponential NoVaS, while the performance of the Limit Model-Free Generalized NoVaS methods are not as good as the Generalized NoVaS under both simulated datasets.

If we look at the results for the three real world financial series in Tables 7–9, the NoVaS still outperforms the benchmark GARCH(1,1) models. Looking at Tables 5 and 6, GARCH(1,1) performs worse when the window size is increasing from 250 to 500, that is, the average coverages are smaller and the average length and standard errors of the predicted intervals are larger. The possible reason might be that the assumption of stationarity might be broken when the window size is too large, like 500. It is also apparent that the performance of Limit Model-Free Generalized NoVaS methods are indistinguishable from that of Generalized NoVaS although the latter appears to be marginally better upon closer comparison. All in all, when the data is non-stationary, the results from Tables 6–9 are supporting the superior performance of time-varying NoVaS in interval prediction against the benchmark time-varying GARCH(1,1). These results are not only consistent with the results of point predictions in Section 2 and Section 3, but also consistent with the findings of the performance of NoVaS without using time-varying processes; see more in Politis (2015), Chapter 10.