2.1. Combining with Normality when the Conditional Variances Are Known

We first consider the case in which the conditional distribution of Y_i given Zⁱ⁻¹ = zⁱ⁻¹ and X_i = x_i is Gaussian for all i ≥ 1 with the conditional variances v_i known. This of course is not realistic in most applications, but it provides a relatively simple case to begin with. To combine the forecasting procedures Δ = {δ₁, δ₂,…} at each time n, we look at their past performances and assign weights accordingly as follows. Let π = {π_j : j ≥ 1} be the prior weights on the procedures, where π_j are positive numbers summing to 1. Let W_j,1 = π_j and for n ≥ 2, let

Note that

for n ≥ 1 and W_j,n depends only on the past forecasts and the corresponding actual realizations of Y. The combined forecasting procedure δ^* is a convex combination of the original procedures:

for n ≥ 1.

Note that in the determination of W_{j, n}, the previous prediction errors are weighted relative to the actual conditional variances v_i. When there is larger uncertainty in Y_i given zⁱ⁻¹ and x_i, prediction error is likely to be larger. The division by v_i accounts for the willingness to tolerate a bigger error in such a case.

Note also that

Thus after each additional observation, the weights on the candidate forecasts are updated. We call the algorithm aggregated forecast through exponential reweighting (AFTER).

This weighting method, as will be seen, leads to a general adaptation risk bound. From (1), the method has a Bayesian interpretation. If we view the weights W_j,n−1, j ≥ 1 as the prior probabilities on the procedures before observing Y_n−1, then W_j,n is the posterior probability of δ_j after Y_n−1 is seen.

For the choice of prior weights {π_j}, if the size of Δ is not large, one can naturally use uniform weights. When the size of Δ is very large, however, one needs to be more subjective, but its effect becomes weaker and weaker as n gets larger.

Related algorithms focusing on worst-case performance have been studied in machine learning (see, e.g., Vovk, 1990; Littlestone and Warmuth, 1994; Cesa-Bianchi, Freund, Haussler, Schapire, and Warmuth, 1997; Haussler, Kivinen, and Warmuth, 1998). Closely related ideas on Bayesian updating, universal coding, and prequential statistical analysis have been considered earlier and since in information theory and statistics (see, e.g., Cover, 1965; Dawid, 1984; Foster, 1991; and recent review articles on universal prediction, Merhav and Feder, 1998; and the minimum description length principle in estimation and prediction, Barron et al., 1998). Results in the machine learning literature focus on bounding the cumulative loss

(or other loss functions) without probabilistic assumptions. Results are usually obtained under the strong condition that the observations Y_i are uniformly bounded. Methods and statistical risk bounds for combining estimation procedures in the context of density estimation, regression, and conditional probability estimation are in Yang (2000a, 2000b, 2000c, 2001) and Catoni (1999).

Condition 0. There exists a constant τ > 0 such that for all i ≥ 1, with probability one we have

This condition means that the individual forecasts are not too far away from the conditional mean relative to the conditional standard deviation. For regression-based forecasts with the regression function bounded below and above by known constants, if the error variance does not diminish, Condition 0 is satisfied when the forecasts are restricted to the same range. Condition 0 is not directly satisfied by the usual time series models such as autoregressive integrated moving average (ARIMA) models. However, like the regression case, for a case where the conditional mean of Y is naturally bounded (e.g., between 0 and 1), restricting the forecasts accordingly can lead to the satisfaction of Condition 0. For instance, for a qualitative threshold autoregressive conditional heteroskedasticity (ARCH) model (Gouriéroux and Monfort, 1992), one may restrict the parameters properly to satisfy Condition 0. It seems quite possible to relax Condition 0 so that specific time series models are allowed without restriction, but we will not pursue that direction in this work. From a practical point of view, in implementing AFTER, we feel that one does not need to worry about this technical condition. Indeed, following the work presented in this paper, Zou and Yang (2003) empirically study the performance of AFTER for ARIMA models and show through simulations and data examples that AFTER has a clear advantage in forecasting compared with popular model selection criteria (used for choosing ARIMA orders) when the instability of model selection cannot be ignored (as is usually the case in most applications).

THEOREM 1. Assume that Condition 0 is satisfied. When the errors are normally distributed, the mean average square risk (relative to the conditional variance) of the combined procedure satisfies

Remarks.

1. Bounds of the type

for combining M procedures based on related methods are given in, e.g., Haussler et al. (1998). But the result is obtained under the very strong assumption that Y_i are uniformly bounded. Such a requirement is usually not appropriate for modeling uncertainty in Y.

2. An upper bound of the type

with a constant C > 1, when taking expectation, is much weaker than the bound given in the theorem.

3. Condition 0 implies that each of the candidate forecasting procedures has scaled mean average square risk upper bounded by τ.

4. Though the risk bound is given under Condition 0, we expect AFTER to perform well more generally. For a good forecasting procedure, (2) is most likely satisfied with high probability, and for a bad procedure, as seen from (1), the procedure's weight drops exponentially fast (if a forecasting procedure is really bad, it will soon be terminated).

From Theorem 1, up to a constant factor and an additive penalty log(1/π_j)/n, the combined procedure achieves the performance of the optimal forecasting procedure among the candidates. In particular, if there are M forecasting procedures then with the uniform prior weight π_j = 1/M, we have

For the nonparametric prediction case,

converges more slowly than order 1/n, and therefore the additive penalty is asymptotically negligible. If additionally when the sample size n increases, the individual forecasts get closer and closer to the conditional means in the sense that

for some positive constants c_i → 0 as i → ∞, then one can show

Thus the performance of the combined forecast is eventually at least half as good as the best individual one (which is unknown). (Note that the condition in (3) is put on all the original procedures and is very restrictive. It is used mainly for technical tractability, and we feel that the upper bound requirement is not really needed for the bad procedures.) The factor of 2 comes from relating the Kullback–Leibler divergence and the Hellinger distance in the technical proof of the theorem.

In Theorem 1, the risks of the forecasts are measured relative to the variability of Y. Under the following conditions (which unfortunately rule out some interesting models such as generalized autoregressive conditional heteroskedasticity (GARCH)), Theorem 1 implies good performance of the combined forecast under the ordinary square loss.

Condition 1. The conditional variances v_n are uniformly upper bounded away from ∞ and lower bounded away from zero in that there exist positive constants A₁ and A₂ such that A₁ ≤ v_n ≤ A₂ a.s. for all n ≥ 1.

Condition 2.

is uniformly bounded above by A₃ with probability one for some constant A₃ > 0.

COROLLARY 1. Assume that Conditions 1 and 2 are satisfied. Under Gaussian errors, the mean average square risk of the combined procedure satisfies

where C is a constant depending only on A₁, A₂, and A₃.

Remark. One does not need to know the constants A₁, A₂, or A₃ to use the AFTER algorithm to combine forecasts.

2.2. Combining with Normality when the Variances Are Estimated

For Theorem 1, the conditional variances v_i are assumed to be known. Similar results can be obtained when the variances are estimated by the individual forecasters. Assume that for each forecasting procedure δ_j, prior to forecasting Y_n, an estimate of v_n, say,

, is obtained based on zⁿ⁻¹ and x_n. For stationary observations a number of variance estimation methods have been proposed for different scenarios. Note that the procedures do not have to use different variance estimators and if some procedures do not provide variance estimates, we can borrow from others.

Condition 3. Assume the variance estimators

are not too far away from the true value: there exist constant 0 < ξ₁ ≤ 1 ≤ ξ₂ < ∞ such that

with probability one for all j ≥ 1 and i ≥ 1.

We modify the AFTER algorithm using the variance estimates. Let W_j,1 = π_j as before and for n ≥ 2, let

Combine the forecasts by

.

THEOREM 2. Assume that the errors are Gaussian and that Conditions 0 and 3 are satisfied. Then the risk of the combined procedure satisfies

where C(ξ₁, ξ₂) = (1/ξ₂ − 1 + log ξ₂)/ξ₁²(1/ξ₂ − 1)².

Regarding the constant C(ξ₁, ξ₂), e.g., when ξ₁ = 1/ξ₂ = ½, C(ξ₁, ξ₂) ≈ 3.1. Note that risks of the variance estimators also appear in the preceding performance bound and thus for this combining method variance estimation is also important. However, we do not need good variance estimation for all the procedures as long as there is an accurate one for a good procedure δ_j. The weighting method and the risk bound suggest that procedures with bad forecasts and/or bad variance estimation receive rather small weights and accordingly have little effect on the performance of good procedures. To have good estimators of v_i, assumptions on the relationship between observations are necessary. For example, one may consider stationary errors that are independent or that have short-range or long-range dependence. Performance of variance estimators can be theoretically justified under reasonable conditions. For various situations, the part of risk in the upper bound due to variance estimation in Theorem 2 may be dominated by the risk of prediction. For example, for prediction based on nonparametric regression with independent or short-range dependent errors, under mild smoothness conditions, the quantity

typically converges at a rate n^−γ for some 0 < γ < 1 but

typically converges at the parametric rate 1/n. In such a case, if the number of individual forecasting procedures in Δ is finite, and if in addition the individual variance estimators converge in the sense that

holds almost surely with ξ_n → 1, then we again have

If no variance estimator is available from the individual forecasters, one can estimate variances based on the previous realizations of Y_i and the corresponding forecasts as follows.

1. Individual variance estimation for the forecasting procedures. Fix a forecasting procedure δ_j. For i = 1, let

be any initial guess. For i > 1, let

Under some conditions (see Proposition 3 in the Appendix), one can show that

converges to zero in probability. Rates of convergence could also be worked out under additional conditions.

2. Estimating variance based on the combined forecasts. The combined forecasting procedure can also be used to estimate variance. For i = 1, let

be any guess. Then use this value as the variance estimate for all the forecasting procedures and accordingly assign weights by (4). For i > 1, let

Other estimators (e.g., based on exponential smoothing) may be more efficient under additional assumptions.

2.3. Combining under Non-Gaussian Errors

When there is a clear indication that the normality assumption on the conditional distributions of Y_n does not hold, the combining weights should be modified. Let h(t), t ∈ R be a probability density function with mean zero and variance 1, which is thought to better capture the error distribution. Non-Gaussian choices of h allow different degrees of heavy-tailedness in the error distributions.

Condition 4. For each pair 0 < s₀ < 1 and T > 0, there exists a constant B = B_s0,T (depending on s₀ and T) such that

for all s₀ ≤ s ≤ s₀⁻¹ and −T ≤ t ≤ T.

Note that ∫ h(y)log(h(y))/(1/s)h((y − t)/s)dy is the Kullback–Leibler divergence between density h and a rescaled and relocated density (1/s) × h((y − t)/s). Under mild conditions, this divergence should behave locally like a squared distance. The condition is satisfied by Gaussian, double-exponential, t (with degrees of freedom bigger than 2), and many other distributions.

We assume that variance estimators are available for the forecasting procedures. Note that the variance estimators given in (5) and (6) in the previous section do not require the normality assumption on the errors; thus they can be used here also. Let

for j ≥ 1 and i ≥ 1. The earlier AFTER algorithm is now modified to compute the weight based on the density function h and the variance estimates. Let W_j,1 = π_j and for n ≥ 2, let

Combine the forecasts by

.

THEOREM 3. Assume that Conditions 0, 3, and 4 are satisfied. Suppose that h used in (7) indeed characterizes the distributions of the errors in the sense that for i ≥ 1, Y_i has the conditional density h((y − m_i)/s_i)/s_i given Zⁱ⁻¹ = zⁱ⁻¹ and X_i = x_i, where

is the conditional standard deviation of Y_i. Then the risk of the combined procedure satisfies

where

.

2.4. Combining under Density Forecasts

Suppose that at each time i, a forecaster δ_j comes up with an estimate, say,

, of p_i, the conditional density of Y_i given Zⁱ⁻¹ = zⁱ⁻¹ and X_i = x_i. The estimate reflects the forecaster's assessment of the uncertainty in Y_i prior to its occurrence. Let

be the mean of the distribution of

. Various methods have been suggested to combine probability distribution estimates (for a review on this topic, see, e.g., Genest and Zidek, 1986). Here our focus is on combining the density estimates to form a composite forecast of Y.

Let W_j,1 = π_j, for n ≥ 2, let

and combine the forecasts by

.

Some technical conditions (to be given subsequently) are needed for part of the next result when

is of the form

, where g_j,i is a (nonrandom) density function with mean zero and variance 1. One hopes that for at least one forecaster δ_j, g_j,i gets closer and closer (in some sense) to p_i as i → ∞.

Let F = {g_θ : θ ∈ Θ} be a class of density functions on the real line. Let f be a given density (f may or may not be in F). Let D(f ∥g) denote the Kullback–Leibler divergence between two densities f and g. We say that f has a projection in F if there exists θ^* ∈ Θ such that for all θ ∈ Θ, D(f ∥g_θ) = D(f ∥g_θ^*) + D(g_θ^*∥g_θ). The density g_θ^* is called the information projection of f in F. When F is an exponential family with

, where φ_i are linearly independent basis functions and ψ(θ) is the normalizing constant, the information projection g_θ^* of f exists and is characterized by ∫ φ_i(x) f (x) dx = ∫ φ_i(x)g_θ^*(x) dx for 1 ≤ i ≤ k (assuming that the integrals are finite) (see Csiszár, 1975).

Now let F_j,i = {g_j,i^{μ, σ} = (1/σ)g_j,i((x − μ)/σ) : −∞ < μ < ∞, σ > 0} be the location-scale family generated by g_j,i for j ≥ 1 and i ≥ 1. We need the following regularity condition for at least one forecasting procedure δ_j.

Condition 5. Suppose that at time i, the true conditional distribution of Y_i has density p_i = h_i((y_i − m_i)/s_i)/s_i, where h_i is a probability density function with mean 0 and variance 1 and is assumed to have an information projection in F_j,i. Assume also that the densities g_j,i, i ≥ 1 satisfy Condition 4 uniformly in the sense that the constants B_s0,T in the condition do not depend on i. Assume further that there exists a constant 0 < c₀ < 1 such that D(g_j,i^{μ, σ}∥g_j,i) ≥ c₀ min(μ² + (σ − 1)²,1) for all i ≥ 1.

As mentioned earlier, the projection requirement in Condition 5 is satisfied if F_j,i are exponential families. The lower bound condition on the Kullback–Leibler divergence within each family is mild and typically satisfied.

THEOREM 4. Assume that Conditions 0 and 3 are satisfied.

Remarks.

Note that because the forms of the conditional densities h_i are unknown, the combining method pays an additional price of order

. If for the forecaster δ_j, g_j,i gets closer and closer to h_i in the Kullback–Leibler divergence, then this penalty approaches zero. It may or may not affect the rate of convergence of the combined procedure, depending on how close g_j,i and h_i are to each other.

In particular, from the result, if there is at least one forecasting procedure that eventually picks h_i with good estimates of the conditional means and standard deviations, the combined procedure based on the AFTER algorithm behaves well. Thus different distribution shapes and assumptions on the dependence relationship between the errors e_i and also on the relationship between Y and the outside information X can be tried out to increase the chance of capturing the true uncertainty in Y. The theorem shows that if one such choice approximates the reality well, the combined forecast will automatically perform well.

2.5. A Consistency Property of the Combining Methods

Recall that the prediction error has a decomposition

. Equivalently, we have

This quantity measures prediction accuracy conditioning on the information available at the time. Note that it is always lower bounded by 1 and should be close to 1 for a good forecast. This motivates the following definition.

DEFINITION 2. A forecasting procedure δ with forecasts

at time i = 1,2,… is said to be consistent if

Note that from the decomposition mentioned previously, δ is consistent if and only if

which roughly means that the forecasts approach the conditional mean m_i (relative to the conditional standard deviation

). This is a desirable property for a forecasting procedure.

PROPOSITION 1. Assume that the conditions for Theorems 1–3 are satisfied. Suppose that there is at least one consistent forecasting procedure in the original list. For the combining methods used in Theorems 2 and 3, we further assume that for at least one consistent forecasting procedure δ_j^*, the variance estimator becomes more and more accurate in the sense that

. Then the combined forecasting procedure is consistent.

The proposition basically indicates that if any original procedure works well, so does the combined procedure in the consistency sense, even if nonconsistent procedures are present. Note that the key characteristic here is that the combining method does not need to know which original procedure is consistent. The conclusion of the proposition follows easily from the results of Theorems 1–3. A similar result can be derived for the combining method in Section 2.4.

2.6. Combining Forecasts without Variance Estimation

The methods in the previous sections require knowledge of the variances or appropriate estimates for determining weights for the candidate forecasts. Sometimes, e.g., when the data are nonstationary, it may be difficult to estimate the conditional variance of Y_n reliably. We present a method in this section that does not require variance estimation. A key idea is borrowed from Catoni (1999).

Condition 6. Assume that the forecasts satisfy that

is bounded above by a positive constant A > 0 with probability one.

Let ψ be a fixed nonnegative convex function with ψ(0) = 0. Here ψ is not necessarily symmetric about 0. When

is used as a loss function, asymmetry of ψ is of interest to reflect the reality that over- and underforecasting might have very different consequences (cf. Newbold and Granger, 1974). Because ψ is convex, for any a₀, there exists a line y = θ_a0(a − a₀) + ψ(a₀) (called a supporting line) such that ψ is above the line. We further assume subsequently that ψ is quadratically off the line and ψ does not change too dramatically in some sense. Let ψ₊′ and ψ₋′ denote the right and left derivatives of ψ, respectively.

Condition 7. There exist constants c > 0 and c > 0 and β > 0 such that for each a₀ ∈ R and a supporting line y = θ_a0(a − a₀) + ψ(a₀) at any a₀, we have

and

for each T > 0.

Condition 7 essentially restricts ψ to be basically quadratic. The slight generality is used so as to include asymmetric function

for some positive constants c₁ ≠ c₂. Such a loss

penalizes prediction error differently when under- or overpredicting Y.

We also need a condition on the conditional distributions of the errors. We do not require knowledge of its shape as long as certain moment generating like functions are uniformly upper bounded.

Condition 8. There exist a constant t₀ > 0 and continuous functions 0 < M₁(t), M₂(t) < ∞ on (−t₀, t₀) such that for all n ≥ 1 and −t₀ ≤ t ≤ t₀, with probability one,

where β > 0 is the same constant as in Condition 7.

Let W_j,1 = π_j and for n ≥ 2, let

where λ is a positive constant to be chosen later. Then the combined forecasting procedure δ^* is defined by

.

THEOREM 5. Assume that Conditions 6–8 are satisfied. Then when λ is chosen small enough, say, 0 < λ ≤ λ₀ for some constant λ₀ depending on c, c, A, M₁, and M₂, the mean average risk (under the loss ψ) of the combined procedure satisfies

In particular, when ψ is the square loss, we have

Note that unlike the previous results, no variance estimation is required here nor is the form of the conditional distribution of Y specified. In addition, the loss (of a constant factor) beyond the additive penalty term is avoided, but under the strong moment generating function conditions. A difficulty with this method is the choice of λ, which can have a dramatic effect on the weights for moderate n. In general, the quantities A, M₁, and M₂ are unknown, and therefore it is hard to determine λ₀. If one is willing to assume a particular form of distribution (e.g., normal), the quantities can be appropriately bounded, and then λ₀ can be determined (for the regression case with independent errors, cf. Catoni, 1999). In contrast, the earlier algorithms do not require such a tuning parameter. Some simulation results in the context of nonparametric regression with independent errors in Yang (2001) suggest that when the true error distribution is double-exponential but weight calculations are based on normal error, estimation accuracy is significantly decreased. On the other hand, when both types of error distributions are considered, the combined procedure behaves as if it knew which error distribution is the correct one. Thus we tend to recommend the use of the earlier methods when one has reasonably strong confidence that some members of the chosen error distribution set closely reflect the true error distribution.

2.7. Potential Applications

2.7.1. Forecasting Based on Regression

Here we use an example to demonstrate how AFTER can be used to obtain forecasts that are adaptive to multiple scenarios.

Suppose that the response variable Y_i and the outside information X_i = (X_i1,…, X_id) have the regression relationship

where m is the unknown regression function assumed to be uniformly bounded and where e_i are Gaussian errors with mean zero and finite variances. Suppose that the explanatory variable X takes values in a bounded cube, which is then taken to be [0,1]^d without loss of generality. In addition, we assume that the errors and X_i are independent of each other and that the unknown design density of X with respect to Lebesgue measure is bounded above and away from zero on [0,1]^d. In practical situations, especially when d is large, it is usually difficult to know the form of m or how the errors are related to each other. Accordingly one may consider different scenarios in the hope that some of them capture the characteristics of the data and result in good forecasts. We consider several models/procedures here. In the first five cases, the errors are assumed to be independent. For such a scenario, as is well-known, prediction of the next response is basically equivalent to estimation of m.

1. Linear regression. Here m is assumed to be of the form

. One can use the least squares method to estimate the parameters. The squared L₂ risk for estimating m is of the parametric rate 1/n if the model holds. One may also consider other plausible parametric models (e.g., including quadratic or higher order and/or cross product terms of the variables).

2. Wavelet regression. For 1 ≤ σ,q ≤ ∞ and α > 0, let B_{q, σ}^{α, d}(C) be the collections of all functions g ∈ L_q [0,1]^d such that the Besov norm satisfies ∥g∥_{Bq, σ^{α, d}} ≤ C (see, e.g., Triebel, 1975; DeVore and Lorentz, 1993). Besov (and similarly defined Triebel's F) classes contain a rich collection of function classes, including Sobolev classes, fractional Sobolev classes, and Hardy classes. The richness of these classes provides a great deal of flexibility (e.g., spatial inhomogeneity) for statistical function estimation. The minimax rate of convergence for estimating a function m ∈ B_{q, σ}^{α, d}(C) under the squared L₂ loss is n^{−2α/(2α+d)} (see, e.g., Donoho and Johnstone, 1998; Yang and Barron, 1999). For d = 1, Donoho and Johnstone (1998) show that wavelet thresholding techniques naturally lead to minimax-rate optimal estimators. Furthermore, the estimators are adaptive in the sense that the hyperparameters σ,q, and α need not be known in advance. For d > 1, one can consider tensor-product wavelet bases and expect similar minimax optimal results.

3. Tensor-product wavelets with different interaction order. Suppose d is large. When the smoothness parameter α is small or moderate, the convergence rate n^{−2α/(2α+d)} for the Besov class B_{q, σ}^{α, d}(C) is slow: the phenomenon known as the “curse of dimensionality.” To overcome the “curse,” one may consider (somewhat enlarged) Besov classes of different interaction orders as follows:

Note that the simplest function class S_{q, σ}^α,1(C) contains additive functions (no interaction). These classes have different effective input dimensions (between 1 and d). By results of Yang and Barron (1999), the minimax rate of convergence under squared L₂ loss for estimating m in S_{q, σ}^{α, r}(C) is n^{−2α/(2α+r)} for 1 ≤ r ≤ d, as suggested by the heuristic dimensionality reduction principle of Stone (1985). Note that when r is small relative to d, the convergence rate is much improved compared with n^{−2α/(2α+d)}. Corresponding to these classes of different interaction orders, one can consider tensor-product wavelets of different interaction orders. With an appropriate thresholding, one expects to obtain an estimator that converges optimally for every S_{q, σ}^{α, r}(C) without knowing the hyperparameters σ, q, α, and r.

4. Neural nets. Let N(C) be the closure in L₂ [0,1]^d of the set of all functions g : R^d → R of the form g(x) = c₀ + [sum ]_i c_i σ(v_i′x + b_i) (where the prime denotes the transpose), with |c₀| + [sum ]_i|c_i| ≤ C, and ∥v_i∥ = 1, where σ is the step function σ(t) = 1 for t ≥ 0, and σ(t) = 0 for t < 0. The minimax rate for estimating m ∈ N(C) under the squared L₂ loss is shown to be bounded between

(see Yang and Barron, 1999). When d is large, the rate is slightly better than n^−1/2 (independent of d), which avoids the “curse of dimensionality.” Estimators at rate O(log n/n^1/2) using finite-dimensional neural network models are in, e.g., Barron (1994).

5. Kernel regression. With kernel methods, universally consistent estimators have been derived under L_q loss without any assumption on the joint distribution of (X,Y) other than the necessary existence of the corresponding moment of Y (see, e.g., Stone, 1977; Devroye and Wagner, 1980). Thus we have a consistent kernel estimator under the squared L₂ loss.

Now suppose that one suspects that the random errors may be correlated. For simplicity, assume that the errors are stationary and normally distributed with Cov(e_i, e_i+j) of order | j|^−γ for some γ > 0. When 0 < γ ≤ 1, the errors are said to be long-range dependent.

6. Prediction under long-range dependence. Suppose that the errors are long-range dependent. Then estimation of the regression function can be severely hampered (see Hall and Hart, 1990; Wang, 1996; Johnstone and Silverman, 1997; Efromovich, 1999; Yang, 1997). However, as long as the dependence parameter is known, prediction is not harder at least in terms of rate of convergence (Yang, 1997). When the dependence is unknown, the prediction problem becomes more complicated. At this point, we are not aware of general results in that direction. We here assume that γ is known and the regression function m has the neural net representation as in case 4. Then one can construct a forecasting procedure converging at rate O(n^−1/2).

Now we combine the forecasts based on the preceding six cases using the method in Section 2.2. If the errors are independent, then when the true regression function is linear in the variables x_i,1 ≤ i ≤ d, the prediction risk of the combined procedure converges at rate 1/n; when m is in B_{q, σ}^{α, d}(C) with α relatively large compared to d, the risk converges at a good rate n^{−2α/(2α+d)}; when m is in S_{q, σ}^{α, r}(C) for some small r, then the risk converges at rate n^{−2α/(2α+r)}; when m is not in any of these cases, but has the neural net representation, then the risk also converges at a good rate O(log n/n^1/2). Even if m is not smooth at all, the combined forecasting procedure is still consistent. If the errors are long-range dependent with the dependence parameter correctly specified in case 6, and if m has the neural net representation, the risk converges at rate O(n^−1/2).

The key point here is that by combining the different forecasts, the mixed forecasting procedure automatically adapts to different scenarios for a good rate of convergence.