2. THE MODEL

Consider the equidistant time series model

where μ is an unknown constant, t_i = i/n, σ(t) > 0 is a smooth, bounded scale (or volatility) function, and {ε_i} is assumed to be a GARCH(r,s) process defined by

(Bollerslev, 1986), where η_i are independent and identically distributed (i.i.d.) N(0,1) random variables, α₀ > 0 and α₁,…,α_r,β₁,…,β_s ≥ 0. Let v(t) = σ²(t) denote the local variance of Y_i. The rescaled time index t_i = i/n is introduced to guarantee that the information increases as n increases and the availability of a consistent estimator of v. Now, model (1) defines indeed a sequence of processes.

Let θ = (α₀,α₁,…,α_r,β₁,…,β_s)′ be the unknown parameter vector. It is assumed that

, which ensures the existence of a unique strictly stationary solution of (2). The practical implementation of a nonparametric estimator

requires the moment condition E(ε_i⁸) < ∞. However, as pointed out by an anonymous referee, the condition of E(ε_i⁴) < ∞ is sufficient for the derivation of the asymptotic results. Necessary and sufficient conditions that guarantee the existence of high-order moments of a GARCH process may be found in Ling and Li (1997), Ling (1999), and Ling and McAleer (2002). It is further assumed var(ε_i) = E(ε_i²) = 1, implying

, to avoid identifiability problems.

The process defined by (1) and (2) is locally stationary in the sense of Dahlhaus (1997), which is a special case of Example 1 given there. Such a model provides a semiparametric extension of the standard GARCH model (Bollerslev, 1986) by introducing the scale function σ(t) into it, where h_i^1/2 stand for the conditional standard deviations of the standardized process ε_i. The total standard deviation at t_i is hence given by σ(t_i)h_i^1/2. For σ(t) ≡ σ₀, model (1) and (2) reduces to the standard GARCH model. Our purpose is to estimate v(t) and h_i separately. If the scale function σ(t) in (1) changes over time, then the assumption of a GARCH model is a misspecification. In this case the estimation of the GARCH model will be inconsistent. It can be shown through simulation that, if a nonconstant scale function is not eliminated, one will obtain

by a fitted GARCH(1,1) model as n → ∞, even when ε_i are i.i.d. Furthermore, in the presence of scale change the estimation of v(t) is also necessary for the prediction. On the other hand, if Y_i follows a GARCH model but model (1) and (2) is used, then the estimation is still

-consistent but with some loss in efficiency due to the estimation of σ(t).

The assumptions of model (1) and (2) can be weakened in different ways. For instance, if the constant mean μ in (1) is replaced by a smooth mean function g, then we obtain the following nonparametric regression with heteroskedastic and dependent errors:

where {ε_i} is a zero mean stationary process. Estimation of the mean function g in model (3) with i.i.d. ε_i is discussed in, e.g., Ruppert and Wand (1994), Fan and Gijbels (1995), and Efromovich (1999). Discussion on the estimation of the scale function in heteroskedastic nonparametric regression may be found in, e.g., Efromovich (1999). This paper focuses on investigating the estimation of σ(t) and θ under model (1) and (2).

3. A SEMIPARAMETRIC ESTIMATION PROCEDURE

Model (1) and (2) can be estimated by a semiparametric procedure combining nonparametric estimation of v(t) and parametric estimation of θ. A linear smoother of the squared residuals will estimate v(t). Let Z_i = (Y_i − μ). Then model (1) can be rewritten as follows:

where X_i = Z_i² and ξ_i = ε_i² − 1 ≥ −1 are zero mean stationary time series errors. Model (4) transfers the estimation of the scale function to a general nonparametric regression problem (for a related idea, see Efromovich, 1999, Sect. 4.3). On the one hand, model (4) is a special case of (3) with g(t) and σ(t) both being replaced by v(t). On the other hand, model (4) also applies to (3) by defining Z_i = Y_i − g(t_i). Hence, the extension of our results to model (3) is expected.

The kernel estimator of conditional variance proposed by Feng and Heiler (1998) will be adapted to estimate v(t). Let y₁,…,y_n, denote the observations. Let

. Let K(u) denote a second-order kernel with compact support [−1,1]. The Nadaraya–Watson estimator of v at t based on

is defined by

where

and b is the bandwidth. And we define

. It is assumed that b → 0, nb → ∞ as n → ∞, which together with other regular conditions ensures the consistency of

. The estimator defined in (5) does not depend on the dependence structure of the errors because

is a linear smoother. It is clear that

if all observations for which |t_i − t| ≤ b are not identical. The bias of

at a boundary point is of a larger order than in the interior because of the asymmetry in the observations. This is the so-called boundary effect of the kernel estimator, which can be overcome by using a local linear estimator (see, e.g., Härdle, Tsybakov, and Yang, 1998). However, as mentioned in Feng and Heiler (1998), a local linear estimator of v may sometimes be nonpositive. Hence, the kernel estimator is more preferable in the current context.

Following Bollerslev (1986), the conditional Gaussian log-likelihood in a parametric GARCH model takes the form (ignoring constants)

The maximizer of L(θ), denoted by

, is not available, because ε_i are unobservable in the current context. Hence we define the approximate log-likelihood by

where

are the standardized residuals given by

The symbols

are used to indicate that, for a given value of θ, h_i(ε;θ) in L(θ) depends on

. Similar to the parameter estimation in the SEMIFAR (semiparametric fractional autoregressive) model (Beran, 1999), θ will be estimated by

, the maximizer of

. Any standard GARCH packet can be used for estimating

. In this paper the S+ GARCH will be used.

obtained in this way is an approximate maximum likelihood estimator (MLE), which may perform differently from

(provided

were available).

5. THE PROPOSED DATA-DRIVEN ALGORITHM

A plug-in bandwidth selector may be developed by replacing the unknowns c_f, I(v²), and I((v′′)²) in (18) with some suitable estimators. At first, it is proposed to estimate c_f by

where

is a nonparametric estimator of E(ε_i⁴). Although explicit formulas of E(ε_i⁴) are known (for common results, see He and Teräsvirta, 1999a; Karanasos, 1999; for results in some special cases, see Bollerslev, 1986; He and Teräsvirta, 1999b), we prefer to use

defined in (21) because the formulas of E(ε_i⁴) are in general too complex. For a GARCH(1,1) model, another simple estimator,

, say, may be defined based on (15) by replacing α₀, α₁, and β₁ with their estimates. Now

perform quite similarly. Assume that a bandwidth b_ε is used for estimating E(ε_i⁴), which satisfies A4 but is not necessarily the same as b. Furthermore, make the following assumption.

A1′. The same as A1 but with E(ε_i⁸) < ∞.

Then the following proposition holds.

PROPOSITION 1. Under Assumptions A1′ and A2–A4 we have

and

where c_f^ε denotes the value of the spectral density of the process ε_i⁴ at the origin.

The proof of Proposition 1 is given in the Appendix.

Remark 1. Equations (22) and (23) show that

is

-consistent, if O(n^−1/2) ≤ b_ε ≤ O(n^−1/4). The optimal bandwidth in a second-order sense, which balances the two terms on the right-hand side of (22), is of order O(n^−1/3). In this paper, we propose to use a bandwidth b_ε = O(n^−1/4) for estimating E(ε_i⁴) so that the estimator is more stable. Note that

is no longer

-consistent if a bandwidth b_ε = O(b_A) = O(n^−1/5) is used. The finally selected bandwidth is not so sensitive to the bandwidth for estimating E(ε_i⁴).

The integral I(v²) can be estimated by

where n₁ and n₂ denote the integer parts of nΔ and n(1 − Δ), respectively, and

is the same as defined in (5) but obtained with another bandwidth b_v, say, that satisfies A4. The following results hold for

.

PROPOSITION 2. Under the assumptions of Proposition 1 we have

and

The proof of Proposition 2 is given in the Appendix.

Remark 2. Note that the dominated orders of the biases and variances of

are the same. Hence similar statements as given in Remark 1 apply for results given in (25) and (26). This is not surprising because both v²(t_i) and ε_i⁴ are related to the fourth moment of the errors.

A well-known estimator of I((v′′)²) is given by

(see, e.g., Ruppert, Sheather, and Wand, 1995), where

is a kernel estimator of v′′ using a fourth-order kernel K₂ for estimating the second derivative (see, e.g., Müller, 1988) and again another bandwidth b_d. Corresponding results as given in Proposition 2 hold for

, for which the following adapted assumptions are required.

A2′. The function v(t) is strictly positive on [0,1] and is at least four times continuously differentiable.

A3′. v′′ is estimated with a symmetric fourth-order kernel for estimating the second derivative with compact support [−1,1].

A4′. The bandwidth b_d satisfies b_d → 0 and nb_d⁵ → ∞ as n → ∞.

PROPOSITION 3. Under Assumptions A1′–A4′ we have

and

The proof of Proposition 3 is omitted because it is well known in nonparametric regression (for results with i.i.d. errors, see, e.g., Ruppert et al., 1995; for results with dependent errors, see, e.g., Beran and Feng, 2002a, 2000b).

Remark 3. The MSE (mean squared error) of

is dominated by the squared bias. The optimal bandwidth for estimating I((v′′)²), which balances the two terms on the right-hand side of (28), is of order O(n^−1/7). With a bandwidth b_d = O(n^−1/7) we have

.

We see that for selecting the bandwidth b we have to choose at first three pilot bandwidths b_ε, b_v, and b_d. This problem will be solved using the iterative plug-in idea (Gasser, Kneip, and Köhler, 1991) with a so-called exponential inflation method (see Beran and Feng, 2002a, 2002b). Let b_j−1 denote the bandwidth for estimating v in the (j − 1)th iteration. Then in the jth iteration, the bandwidths b_ε,j = b_v,j = b_j−1^5/4 and b_d,j = b_j−1^5/7 will be used for estimating E(ε⁴), I(v²), and I((v′′)²), respectively. These inflation methods are chosen so that b_ε,j and b_v,j are both of order O_p(n^−1/4) and b_d,j is of the optimal order O_p(n^−1/7), when b_j−1 is of the optimal order O_p(n^−1/5). By an iterative plug-in algorithm the unknown constants in the pilot bandwidths can be simply omitted. Furthermore, we also need to choose a starting bandwidth b₀. In the current context, b₀ should satisfy A4 because we have to estimate θ in the first iteration. Theoretically, a bandwidth b₀ = O(n^−1/5) is more preferable. Our experience shows that b₀ = 0.5n^−1/5 is a good choice. Detailed discussions on this topic may be found in the next two sections, especially in Section 6.3.

The proposed data-driven algorithm processes as follows:

1. Start with the bandwidth b₀ = c₀ n^−1/5 with, e.g., c₀ = 0.5.

2. In the jth iteration

3. Increase j by one and repeatedly carry out step 2 until convergence is reached or until a given maximal number of iterations has been completed. Put

.

The condition |b_j − b_j−1| < 1/n is used as a convergence criterion of

, because such a difference is negligible. The maximal number of iterations is put to be 20. In this algorithm,

is estimated using b_j−1 as for

because we do not have a proper bandwidth selector for estimating θ. The asymptotic performance of

is quantified by the following theorem

THEOREM 4. Assume that A3 and A1′–A3′ hold and that I((v′′)²) ≠ 0. Then we have

The proof of Theorem 4 is given in the Appendix. Note that A4 and A4′ are automatically satisfied. The second O(n^−2/5) term on the right-hand side of (31) is due to the error in

caused by the bias in

, which is indeed negligible compared with the first term.

The proposed algorithm is coded in an S-Plus function called SEMIGARCH. A practical restriction 1/n ≤ b ≤ 0.5 − 1/n is used in the program for simplicity. Four commonly used kernels, namely, the uniform, the Epanechnikov, the bisquare, and the triweight kernels (see, e.g., Müller, 1988), are built into the program. As a standard version we propose the use of the Epanechnikov kernel with Δ = 0.05 and c₀ = 0.5, which will be used in the next two sections.

Remark 4. Note that b_A is not well defined if I((v′′)²) = 0 implying v′′(t) ≡ 0. However, the SEMIGARCH model also applies to this case. In particular, the proposed algorithm does work if Y_i follow a GARCH model. Now it can be shown that, theoretically, b_j → O_p(1) as j → ∞. Following the context after Theorem 3,

has the same asymptotic properties as by a GARCH model because

. And

is

-consistent with some loss in the efficiency compared with a parametric estimator, provided that no maximal number of iterations is given, because (nb_j)⁻¹ → O_p(n⁻¹) now.

6. THE SIMULATION STUDY

6.1. Design of the Simulation

To show the practical performance of our proposal, a simulation study was carried out. In the simulation study, ε_i were generated using the simulate.garch function in S+GARCH following one of the two GARCH(1,1) models.

Model 1 (M1). ε_i = η_i h_i^1/2, h_i = 0.6 + 0.2ε_i−1² + 0.2h_i−1 and

Model 2 (M2). ε_i = η_i h_i^1/2, h_i = 0.15 + 0.1ε_i−1² + 0.75h_i−1.

The y_i are generated following model (1) with μ ≡ 0 and one of the three scale functions:

The terms v₁(t) and v₂(t) are quite similar, and they are designed following the estimated scale function in the daily DAX 100 returns. The scale change with v₂ is stronger than that with v₁. It is strongest with v₃. To this end see the bandwidths required for estimating them given in Table 2. The scale function σ₂(t) may be found in Figure 2b, which follows. To confirm the statements in Remark 4, a constant scale function v₀(t) = σ₀²(t) ≡ 16 is also used. The simulation was carried out for three sample sizes n = 1,000, 2,000, 4,000. For each case 400 replications were done. For each replication, three GARCH (1,1) models were fitted to ε_i, the data-driven

. The estimators of α₁ and β₁ are denoted by

, respectively. For v₀ we have

. Here,

are used as a benchmark. Note in particular that the estimated parameters may sometimes be negative using the S+GARCH.

6.2. Results of the Simulation Study

To give a summary of the performance of

, and to compare them with

, the empirical efficiency (EFF) of an estimator w.r.t. the corresponding one estimated from ε_i is calculated. For instance,

These results are listed in Table 1. The difference between two related EFFs, e.g.,

, in a given case may be thought of as the gain by using the SEMIGARCH model. Table 1 shows that the EFFs of

seem to tend to 100%, whereas those of

seem to tend to zero, as n → ∞. Hence, the gains seem to tend to 100% as n → ∞. The EFFs of

under M2 are relatively low. In particular, for n = 1,000, the EFFs of

in the two cases of M2 with v₁ and v₃ are even smaller than those of

; i.e., the gain in these cases is slightly negative. This shows that n = 1,000 is sometimes not large enough for estimating the scale function when β₁ is large.

Empirical efficiencies (%) of the estimated parameters

Box plots of the 400 replications of

for n = 1,000 are shown in Figures 1a–1f, where the symbols E1, E2, and E3 denote estimators obtained from

, respectively. Those for n = 2,000 and n = 4,000 are omitted to save space. The simulation results show that

perform in general quite well. One clear problem arises with

under M2 with n = 1,000. Now, both the variance and the bias of

are strongly affected by some extremely small estimates (see Figures 1m–1p). This is due to the nonrobustness of the bandwidth selection. Hence, it is worthwhile to develop a robust procedure to improve the poor performance of

for small n. The quality of

is clearly improved as n increases. In particular, the estimation becomes more and more stable. Detailed statistics (in the first version of this paper) show that the standard deviations of

seem to converge at the same rate as those of

, but their biases converge a little more slowly. This confirms the results of Theorem 3. The simulation results show clearly that, in the case with scale change,

are inconsistent as a result of their biases. The situation become worse as n increases. In particular, we can see that

will tend to one as n → ∞, no matter how large β₁⁰ is. However, if there is no scale change the estimators

should of course be used. It is hence helpful to test whether or not the estimated scale function is significant. For the data examples given in the next section it is proposed to carry out such a test based on simulation.

Box plots of (E3), respectively, with n = 1,000, where the horizontal lines show the true values.

Now let us consider the quality of

. The sample means, standard deviations, and square roots of the MSEs of

together with the true asymptotic optimal bandwidths b_A are given in Table 2. Note that b_A and the MSE in cases with v₀ are not defined. Kernel density estimates of

(omitted to save space) show that the performance of

is satisfactory. In all cases the variance of

decreases as n increases. It is also true for the bias in most of the cases. Both the variance and the bias of

depend on the scale function and the model of the errors. For two related cases, the variance of

under M1 is smaller than that under M2. Generally, the stronger the scale change, the larger the variance of

. The bias of

by v₁ is always negative, and it is always positive by v₃. The bandwidth for v₂ is easiest to choose. The choice of the bandwidth by v₃ is in general easier than that by v₁, except for the case of M2 with n = 1,000. In this case, the detailed structure of v₃ may sometimes be smoothed away because of the large variation caused by the GARCH model. This shows again that n = 1,000 is sometimes not large enough for distinguishing the CH and the scale change.

Statistics on the selected bandwidth

Remark 5. As suggested by a referee, the performance of the proposed procedure for the cases with highly persistent GARCH effect is investigated through an additional simulation under a third model, model 3 (M3), with α₁ = 0.07 and β₁ = 0.87 and without trend. As expected, the proposed procedure does not work well for n = 1,000 because the variance of

and in particular that of

are too large as a result of some extreme estimates. This shows again that a robust estimation procedure should be developed. For n ≥ 2,000, the procedure works well. The empirical efficiencies are a little lower than those for M2. Detailed results of this additional simulation are omitted to save space.

Remark 6. In this paper, the bandwidth is selected by minimizing the dominant part of the MISE of

. In a semiparametric context, the performance of the bandwidth selection and the resulting parameter estimation may be improved if a plug-in algorithm that takes the MSE of

into account is developed. For this purpose a more detailed formula of the MSE of

is required, and one has to develop a suitable procedure to estimate the MSE. This is still an important open question and will be discussed elsewhere.

6.3. Detailed Analysis of Two Simulated Examples

In the following discussion, two simulated data sets are selected to show some details. The first example (called Sim 1) is a typical example of the replications under M2 with the scale function σ₂(t) and n = 2,000. The observations y_i, i = 1,…,2,000, are shown in Figure 2a. For Sim 1 we have

by starting with any bandwidth 3/n ≤ b₀ ≤ 0.5 − 1/n; i.e.,

does not depend on b₀ if b₀ is not too small. The σ₂(t) (solid line) and

(dashed line) are shown in Figure 2b. Figure 2c shows the standardized residuals

, which look stationary. The estimated GARCH(1,1) models are

Estimation results for the first simulated data set.

for y_i and

for

. For model (32) we have

, so that the fourth moment of this model does not exist. On the opposite model (33) has finite moments until at least twelfth order as for the underlying GARCH model. The estimated SEMIGARCH conditional and total standard deviations, i.e.,

, are shown in Figures 2d and 2e. The true conditional and total standard deviations of y_i, i.e., (h_i)^1/2 and σ₂(t_i)(h_i)^1/2, are shown in Figures 2f and 2g. Figure 2h shows the estimated GARCH conditional (in this case also the total) standard deviations (h_i^y)^1/2. The analysis of Sim 1 shows the following results.

(1) If a standard GARCH model is used, the scale change will be wrongly estimated as a part of the CH. Furthermore, the total variance tends to be overestimated when it is large and underestimated when it is small (compare Figures 2g and 2h). This phenomenon is mainly due to the overestimation of

and will be called the (volatility) inflation effect of the GARCH model in the presence of scale change.

(2) Following the SEMIGARCH model, both the conditional heteroskedasticity and the scale change are well estimated. The estimated SEMIGARCH total variances are quite close to the true values and are more stable and accurate than those following the standard GARCH model (compare Figures 2e and 2h). The errors in

are caused by the errors in these two estimates, and both of them can be clearly reduced if more dense observations are available, e.g., by analyzing high-frequency financial data. The MSE of the estimated total variances are 0.687 for the SEMIGARCH and 4.979 for the standard GARCH models; the latter is about seven times as large as the former.

Furthermore, (h_i^y)^1/2 shown in Figure 2h (see also Figure 3f) exhibit a clear signal of covariance nonstationarity, a property not shared by the true and the estimated SEMIGARCH conditional standard deviations.

The estimation results for the S&P 500 returns.

The second simulated data set (called Sim 2) is one of the replications under M1 with v₃ and n = 1,000, which is chosen to show that sometimes the selected bandwidth will be wrong if b₀ is too small or too large. That is, a moderate b₀ should be used as proposed in Section 5. For this data set we have either

if b₀ < 0.020. On the other hand, we have

, the largest allowed bandwidth in the program, if b₀ > 0.262. For any starting bandwidth b₀ ∈ [0.021, 0.262] a bandwidth

will be selected. Now,

does not depend on b₀. Note that the proposed default starting bandwidth b₀ = 0.5n^−1/5 = 0.126 lies in the middle part of the interval [0.021, 0.262]. In case when it is doubtful, if the selected bandwidth with b₀ = 0.5n^−1/5 is the optimal one, we recommend that the user try with some different b₀'s and choose the most reasonable

from all possible selected bandwidths by means of further analysis (see Feng, 2002).

7. APPLICATIONS

In this section the proposal will be applied to the log-returns of the daily S&P 500 and DAX 100 financial indexes from January 3, 1994, to August 23, 2000. For the S&P 500 returns shown in Figure 3a we have

(for any b₀ ≥ 0.075). The fitted GARCH models are

for y_i and

for

. As before, for model (34) we have

so that the fourth moment of this model does not exist. Model (35) has finite moments until twelfth order. To test whether the estimated trend is significantly nonconstant, 400 replications were generated following model (35) with the corresponding sample variance and without trend. The scale function was then estimated with the bandwidth b = 0.183 from each replication. Symmetric Monte Carlo confidence bounds that covered 95% or 99% of all estimated trends were calculated and are shown in Figure 3b together with the sample standard deviation (0.0099) and the estimated scale function

. We see that there is significant scale change in this data set. Furthermore, both

in model (35) are strongly significant. That is, this series has simultaneously significant scale change and CH. Figures 3c–3f show

, the SEMIGARCH conditional standard deviations

, the SEMIGARCH total standard deviations

, and the GARCH conditional standard deviations (h_i^y)^1/2. Comparing Figures 3e and 3f we see again that the estimated total variances following the SEMIGARCH model are more stable and those following the GARCH model are inflated.

For the DAX 100 returns we have

(for any b₀ ≥ 0.075). The fitted GARCH models are

for y_i and

for

. The condition for the existence of the fourth moment of model (36) is slightly satisfied, but the eighth moment of this model does not exist. Again, model (37) has finite moments until twelfth order. The S&P 500 and DAX 100 returns series perform quite similarly, and the conclusions on the former given previously apply to the latter.

Now, we will compare the performance of the GARCH and SEMIGARCH by predicting future volatility. The GARCH unconditional variance,

, say, is calculated following (34) or (36). For the SEMIGARCH,

is used as the unconditional variance in the near future. The predicted (expected) conditional standard deviations

following the GARCH and

following the SEMIGARCH, k = 1,2,…,100, for the S&P 500 and DAX 100 returns are shown in Figure 4 together with

. Note that, the conditional standard deviations by both series at the right end are lower than

. Consequently,

increase for both series. The

look quite reasonable and converge to

quickly. However,

in both cases seem to be underestimated, because of the inflation effect mentioned previously. Furthermore,

converge very slowly to some wrongly estimated limits. The sample standard deviation for the S&P 500 returns is 0.0099. Following (34) we have

, which is clearly overestimated as a result of the instability of this model. For the DAX 100 returns,

is about equal to its sample value, which is, however, clearly lower than the locally unconditional standard deviation at t = 1. There are two problems if the fitted parametric GARCH models from these data sets are used for predicting future volatility: (1) the unconditional variance at the current end was wrongly estimated; and (2) the predicted conditional variance converges very slowly, because these models only have finite moments of low orders. Both of these problems were overcome by applying the SEMIGARCH model.

Predicted standard deviations (middle dashes) and (solid line) together with their limits following the GARCH (short dashes) and following the SEMIGARCH (long dashes) for (a) the S&P 500 and (b) the DAX 100 returns.