Autoregressive time series prediction by means of fuzzy inference systems using nonparametric residual variance estimation

Autoregressive Time Series Prediction by Means of Fuzzy Inference Systems Using Nonparametric Residual Variance Estimation Federico Montesino Pouzols ∗, Amaury Lendasse Department of Information and Computer Science, Helsinki University of Technology P.O. Box 5400 FI-02015 HUT Espoo, Finland Angel Barriga Barros bDepartment of Electronics and Electromagnetism University of Seville Avda. Reina Mercedes s/n, E-41012 Seville, Spain. ∗Corresponding author. Phone: +358-9-451-5237. FAX: +358-9-451-3277. Email addresses: [email protected] (Federico Montesino Pouzols), [email protected] (Amaury Lendasse), [email protected] (Angel Barriga Barros). Preprint submitted to Elsevier 3 November 2009 Abstract We propose an automatic methodology framework for short- and long-term prediction of time series by means of fuzzy inference systems. In this methodology, fuzzy techniques and statistical techniques for nonparametric residual variance estimation are combined in order to build autoregressive predictive models implemented as fuzzy inference systems. Nonparametric residual variance estimation plays a key role in driving the identification and learning procedures. Concrete criteria and procedures within the proposed methodology framework are applied to a number of time series prediction problems. The learn from examples method introduced by Wang and Mendel (W&M) is used for identification. The Levenberg-Marquardt (L-M) optimization method is then applied for tuning. The W&M method produces compact and potentially accurate inference systems when applied after a proper variable selection stage. The L-M method yields the best compromise between accuracy and interpretability of results, among a set of alternatives. Delta test based residual variance estimations are used in order to select the best subset of inputs to the fuzzy inference systems as well as the number of linguistic labels for the inputs. Experiments on a diverse set of time series prediction benchmarks are compared against least-squares support vector machines (LS-SVM), optimally-pruned extreme learning machine (OP-ELM), and k-NN based autoregressors. The advantages of the proposed methodology are shown in terms of linguistic interpretability, generalization capability and computational cost. Furthermore, fuzzy models are shown to be consistently more accurate for prediction in the case of time series coming from real-world applications Key words Fuzzy inference systems, Time series prediction, Nonparametric regression, Supervised learning, Nonparametric residual variance estimation 2 1 Introduction Time series prediction and analysis in general is a recurrent problem in virtually all areas of natural and social sciences as well as in engineering. In the time series prediction field, prediction accuracy is not the only major goal. Understanding the behavior of time series and gaining insight into their underlying dynamics is a highly desired capability of time series prediction methods [52]. In the past, conventional statistical techniques such as AR and ARMA models have been extensively used for forecasting [6]. However, these techniques have limited capabilities for modeling time series data, and more advanced nonlinear methods including artificial neural networks have been frequently applied with success [9]. Fuzzy logic based modeling techniques are appealing because of their interpretability and potential to address a broad spectrum of problems. In particular, fuzzy inference systems exhibit a combined description and prediction capability as a consequence of their rule-based structure [50]. The application of fuzzy inference systems to time series modeling and prediction dates back to [51], in which the authors develop the well known learn from examples identification algorithm for fuzzy inference systems and use the Mackey-Glass time series as a validation case. Nevertheless, despite its good performance in terms of accuracy and interpretability, fuzzy inference systems have seen little application in the field of time series prediction as compared to other nonlinear modeling techniques such as artificial neural networks and support vector machines. The methodology proposed in this paper is intended to apply to crisp time series, i.e. those time series consisting of crisp values, as opposed to other kinds of values, such as interval and fuzzy values. That is, we propose here an automatic methodology framework to perform autoregressive prediction of crisp time series by means of fuzzy 3 inference systems using nonparametric residual variance estimation [37]. We will call fuzzy autoregressors those autoregressors implemented as fuzzy inference systems. This is not to be confused with what is usually called fuzzy regression in the literature [8]. When developing fuzzy inference systems for time series prediction, many questions remain still open: how many and what inputs to the inference system must be defined? To what extent the theoretical universal approximation capability of fuzzy systems is achieved with existing techniques? What are the best fuzzy techniques for these tasks? How to perform long-term prediction? In practice, one finds two problems when building a fuzzy model for a time series: choosing variables or inputs to the inference system, and identifying the structure of the system (linguistic labels and rule base). Once these steps have been accomplished, the fuzzy model can be tuned through supervised learning techniques. We propose an automatic methodology framework to address these two problems using fuzzy techniques and nonparametric residual variance estimation techniques in an intertwined manner. The first problem can be addressed by means of a priori feature selection techniques based on nonparametric residual variance estimation, which also provide an estimate of the error of the most accurate nonlinear model that can be built without overfitting. The second problem is addressed by data-driven techniques for identification of fuzzy systems from numerical examples [16], such as the algorithm by Wang and Mendel (W&M) [50,51] and the fuzzy identification algorithm based on clustering by Chiu [12,36], two well established methods among the many alternatives proposed throughout the years [23,19,43,34,2,26,44]. This paper also addresses a recent challenge in the field of time series prediction: long-term prediction (as a generalization to short-term prediction), for which lack of 4 information and accumulated errors pose additional difficulties. Furthermore, time series coming from real world applications, in addition to series generated either numerically or under controlled laboratory conditions are analyzed. With these premises we propose a methodology for building simple, interpretable yet highly accurate fuzzy inference models. These will be compared against least-squares support vector machines (LS-SVM) [48], a well established method in the field of time series prediction, that has been shown to be highly accurate. For further comparison, we will also show the results obtained using optimally-pruned extreme learning machine (OP-ELM) and k-nearest neighbors (k-NN) models. The article is organized as follows. The next section outlines a nonparametric residual variance estimation method that will be used for both variable and proper model complexity selection. In section 3we propose a methodology framework and one concrete implementation that uses well known algorithms for identifying and optimizing fuzzy inference systems. Section 4illustrates the methodology through a case study. Finally, sections 5and 6present and further discuss experimental results for a number of time series benchmarks from diverse fields of application. In appendices Aand Bwe provide further comparisons of different modeling techniques. 2 Nonparametric Residual Variance Estimation: Delta Test Nonparametric residual variance estimation (or nonparametric noise estimation, NNE) is a well-known technique in statistics and machine learning, finding many applications in nonlinear modeling [24]. NNE methods can be applied to recurrent problems such as variable and model structure selection. These methods are not however in widespread use in the machine learning community as most work has been done to date within the statistics community. 5 Delta Test (DT), introduced for time series in 1994 [41], is a NNE method, i.e., it estimates the lowest mean square error (MSE) that can be achieved by a model without overfitting the training set [24]. Given Nmultiple input-single output pairs, (¯xi, yi)∈RM×R, the theory behind the DT method considers that the mapping between ¯xiand yiis given by the following expression: yi=f(¯xi) + ri, where fis an unknown perfect fitting model and riis the noise. DT is based on hypothesis coming from the continuity of the regression function. When two inputs x and x′are close, the continuity of the regression function implies that the corresponding outputs, f(x)and f(x′)will be close enough. When this implication does not hold, it is due to the influence of the noise. Let us denote the first nearest neighbor of the point ¯xiin the set {¯x1,...,¯xN}by ¯xNN . Then the DT, δ, is defined as follows: δ=1 2N N X i=1 yNN(i)−yi 2, where yNN(i)is the output corresponding to ¯xNN(i). For a proof of convergence, refer to [28,29]. DT is an unbiased and asymptotically perfect estimator with a relatively fast convergence [29] and is useful for evaluating nonlinear correlations between two random variables, namely, input-output pairs. DT can be seen as part of a more general NNE framework known as the Gamma Test [24]. Despite the simplicity of DT, it has been shown to be a robust method in real world applications [28]. This method will be used in the next sections for a priori input selection. 6 3 Methodology Framework for Time Series Prediction with Fuzzy Inference Systems Consider a discrete time series as a vector, ¯y=y1, y2,...,yt−1, yt, that represents an ordered set of values, where tis the number of values in the series. The problem of predicting one future value, yt+1, using an autoregressive model (autoregressor) with no exogenous inputs can be stated as follows: ˆyt+1 =fr(yt, yt−1,...,yt−M+1), where ˆyt+1 is the prediction of model frand Mis the number of inputs to the regressor, i.e., the regressor size. Predicting the first unknown value requires building a model, fr, that maps regressor inputs (known values) into regressor outputs (predictions). When a prediction horizon higher than 1 is considered, the unknown values can be predicted following two main strategies: recursive and direct prediction. The recursive strategy applies the same model recursively, using predictions as known data to predict the next unknown values. For instance, the third unknown value is predicted as follows: ˆyt+3 =fr(ˆyt+2,ˆyt+1, yt, yt−1,...,yt−M+3). Recursive prediction is the most simple and intuitive strategy and does not require any additional modeling after an autoregressor for 1 step ahead prediction is built. However, recursive prediction suffers from accumulation of errors. The longer the prediction horizon is, the more predictions are used as inputs. In particular, for prediction horizons greater than the regressor size, all inputs to the model are predictions. Direct prediction requires that the process of building an autoregressor be applied for 7 each unknown future value. Thus, for a maximum prediction horizon H,Hdirect models are built, one for each prediction horizon h: ˆyt+h=fh(yt, yt−1,...,yt−M+1),with 1≤h≤H. While building a prediction system through direct prediction is more computationally intensive (as many times as values are to be predicted) it is also straightforward to parallelize. As opposed to recursive prediction, direct prediction does not suffer from accumulation of prediction errors. In this paper, we follow the direct prediction strategy. In order to build each autoregressor, a fuzzy inference system is defined as a mapping between a vector of crisp inputs and a crisp output. Let us rename the inputs yt, yt−1,...,yt−M+1 as y1,...,yMfor simplicity. This way, assuming all (M) inputs are used, the fuzzy autoregressor for prediction horizon hcan be expressed as a set of Nhfuzzy rules of the following form: Rh i:IF y1is Li,h 1AND y2is Li,h 2AND ... AND yMis Li,h MTHEN ˆyt+h←µRh i, where i= 1,...,Nh, and the fuzzy sets Li,h j∈ {Lh j,k}, k = 1 ...,nh j, j = 1,...,M,with nh jbeing the number of linguistic labels defined for the jth input variable. Li,h jare the fuzzy sets representing the linguistic terms used for the jth input in the ith rule of the fuzzy model for prediction horizon h.µRh iare the consequents of the rules and can take different forms. For example, in a system with two inputs, if Li,h 1is renamed LOW1and Li,h 2is renamed HIGH2, the ith rule for horizon 1, R1 i, would have the following form: IF ytwas LOW1AND yt−1was HIGH2THEN ˆyt+h←µRh i. Depending on the fuzzy operators, inference model and type of membership functions (MFs) employed, the mapping between inputs and outputs can have different formulations. In principle, the methodology proposed in this paper can be applied for 8 any combination of types of MFs, operators and inference model, but the selection can have a significant impact on practical results. As a concrete implementation for this paper, we use the minimum as T-norm for conjunction operations, Gaussian MFs for inputs, singleton outputs, and product inference of rules. Defuzzification is performed using the fuzzy mean method, i.e., zero-order Takagi-Sugeno systems are defined. Thus, the result of the inference process is a weighted average of the singleton consequents. This inference scheme was chosen in order to keep systems as simple and interpretable as possible. Therefore, in this particular case a fuzzy autoregressor for prediction horizon hcan be formulated as follows: Fh(¯y) = Nh X i=1 µRh i·min 1≤j≤MµLi,h j(yj) Nh X i=1 min 1≤j≤MµLi,h j(yj) ,(1) where Nhis the number of rules in the rule base for horizon h,µRh iare singleton output values, and µLi,h jare Gaussian MFs for the inputs. Thus, each fuzzy set defined for the input linguistic terms, Lh j,k (for horizon h, and the kth term defined for the jth input), is characterized by an MF having the following form: µLh j,k =exp h−(yj−cj,k,h)2/2σ2 j,k,hi, k = 1 ...,nh j, j = 1, . . . , M, h = 1,...,H, where cj,k,h are the centers of the MFs and σj,k are the widths. Together with the number of rules of a system, the total number of MFs can be seen as a measure of the complexity of a fuzzy inference system, or the structure of the equivalent artificial neural network. If the same nh jnumber of linguistic terms is set for every input, then the total number of input MFs per prediction horizon is M·nh j. The problem of building a regressor can be precisely stated as that of defining a proper 9 4.1 Variable Selection Following the flow depicted in figure 1, in the first stage of our methodology, DT is performed on the training set for all the possible variable selections (210 −1) and the one that leads to the lowest DT estimate is chosen. This process is performed independently for each prediction horizon. The number of selected variables is shown in figure 3. Between 3 and 5 variables are selected out of a maximum of 10. Thus, the employment of DT based variable selection leads to a significant decrease of the complexity of the fuzzy inference systems in terms of number of inputs. This fact, in turn, relieves the curse of dimensionality problem. Thus, for a given maximum regressor size, an initial input selection stage allows for a better fitting of the model. We should note that a maximum regressor size larger than 10 could be considered. When there are enough data samples, a larger regressor size could be expected to provide accuracy improvements. However, for sizes above 15 or 20 approximately, an exhaustive search becomes too computationally expensive and finally unfeasible with current computational resources. A regressor size of 10 has thus been selected as a twofold heuristic compromise. First, the whole space of possible selections can be explored within a reasonable amount of time (approximately 1 hour for 50 models in a current general purpose computer). Second, after variable selection, the number of inputs is sufficiently small so that the curse of dimensionality problem in nonlinear models does not have a severe impact. Larger regressor sizes, for which the DT estimate is lower, usually lead however to little improvement or even poorer performance of the models. We thus, have selected 10 as an initial regressor size. As we will see later on in this section, in the case of fuzzy inference models this leads to systems with a number of inputs and rules sufficiently small so as to be easy to read by humans. The effect that different maximum regressor sizes can have on model performance will be illustrated 16 with several cases in section 5. 4.2 Identification and Tuning As second stage, once input variables have been selected, an iterative identification and tuning process is carried out in three substages, as shown in figure 1. In the first substage, the W&M algorithm is applied to the training set in order to identify fuzzy inference systems. These models are then tuned in the second substage through supervised learning using the Levenberg-Marquardt algorithm over the training set. The process is repeated for increasing numbers of linguistic labels (or MFs) per input, starting from 2. Within this iterative process, in the third substage (complexity selection) the DT estimate is used to check whether the best possible approximation has been achieved, i.e., the right compromise between model complexity and training error has been found. For the horizon 1 regressor, table 1shows the number of rules identified for different numbers of linguistic labels per input (between 2 and 15). Training and validation errors are shown as well. The two columns labeled “before tuning” show the errors for the fuzzy systems as identified by means of the W&M algorithm, while the columns labeled “after tuning” show the errors for the systems tuned by means of supervised learning. After the tuning substage, there is a considerable accuracy improvement. In particular, it can be seen that tuned systems with a low number of rules perform better than untuned systems with a much greater complexity. Thus, the supervised learning substage also contributes to reducing model complexity. We also note that systems with a low number of linguistic labels per input (particularly between 2 and 5) are only rough approximators before tuning. However, after the tuning substage their accuracy is improved significantly while keeping the same rule base. This 17 fact suggests that the rule bases correctly reflect the underlying dynamics of the series, though tuning the membership function parameters is no doubt required in order to build accurate models with such a low number of linguistic labels. Within the methodology proposed here, in this case the identification and tuning stage proceeds as follows. First, in the identification substage a fuzzy inference system is identified using the W&M algorithm with 2 MFs per input. The system consists of 6 rules. After the tuning substage, the system yields a training MSE higher than the DT estimate. Thus, this system is rejected in the complexity selection substage. Then, the process is repeated for a system with 3 MFs per input, for which 15 rules are identified. In this case, the system yields an MSE lower than the DT estimate after being tuned. Thus, it qualifies as a proper model in the complexity selection substage and the identification and tuning stage finishes 2. 4.3 Interpretability Issues and Examples The number of MFs or linguistic labels defined for each input has significant consequences on accuracy and interpretability. If it is too low, the accuracy of the system will not suffice. If it is high, the linguistic labels will become too specific, and the number of fuzzy rules identified can be overwhelming, since it can grow up to the product of the number of MFs of every input. Let us now consider the interpretability issues that arise after the tuning substage. Fuzzy inference systems are inherently comprehensible, specially when the rules are defined by human experts. However, when rules are automatically identified from data and optimization methods are applied, as is the case of the proposed methodology, 2In practice, as detailed later on in this section, we use a certain error threshold based on the DT estimate rather than the estimate itself. 18 interpretability cannot be guaranteed in general [7]. Since the identification substage is implemented using the W&M method, the initial input MFs, of Gaussian type, define a grid uniform partition of the input domain. The output MFs, of singleton type. correspond to output values identified on the training data. The W&M method also guarantees the consistency of the rule base. In a system of this kind, linguistic labels can be easily assigned to each input MF by domain experts, or simply by using the common “LOW”, “MEDIUM”, etc. labels. The meaning of singleton output values is evident as well. In addition, the rules identified are of the if-then type with only conjunction operations in the premise. Thus, these systems can be expected to be easy to read and potentially lead to a physical interpretation. However, the tuning substage consists in changing the parameters of the input and output MFs with the objective of finding the lowest MSE. After the tuning stage, interpretability could be thus severely compromised. A variety of approaches and methods have been proposed to improve or guarantee interpretability to some extent [7]. Nonetheless, in our application it can be found in practice that the parameter changes do not modify the initial uniform partition of the input space to an extent significant for approximate, linguistic human interpretation. On the one hand, changes in the output MFs should not decrease interpretability in an automated methodology. On the other hand, the changes in the input MFs do not modify the initial partition in a severe manner. The shape and distribution of the tuned MFs of the regressor for horizon 1 are shown in figure 4. In this case, the widths change by 12.9% on average, while the centers of the MFs are shifted on average by 3.1% of the range of the series, with respect to the initial grid uniform partition. More general results for the set of series analyzed in this paper are given in appendix B, table 6. For the 1 step ahead regressor, considering the notation for discrete time series 19 introduced in section 3, three input variables are selected to predict yt+1:yt,yt−2and yt−7. This selection indicates that the next weekly temperature depends only on the values of the last week as well as 2 and 7 weeks before. In addition, the relations between the inputs and the output can be interpreted linguistically. Three linguistic terms are defined for each input, as shown in figure 4, that can be thought as “LOW”, “MEDIUM” and “HIGH”, represented by Gaussian MFs in the range of observed values, [18.9,29.2]. For instance, a rule that has the highest output temperature as consequent reads as follows: IF yt−7was HIGH1AND yt−2was HIGH2AND ytwas HIGH3THEN yt+1 ←“29.2°C”, where “29.2°C” is used as linguistic label for a singleton output centered at 29.2. Of course, this example can be regarded only as a simple particular case. The procedure required to provide a physical interpretation of the models is case dependent to a great extent. In the interpretation process, additional techniques for simplification of fuzzy inference systems can be of considerable help [4]. For instance, if the system is pruned in order to keep only the six best rules with respect to the training set, a system with a test MSE only 6.3% higher than that of the original system is obtained. Let us now illustrate with an example the process of fuzzification, inference and defuzzification. Consider the application of the previous fuzzy rule to a data point in the test series such that yt+1 = 28.8has to be predicted from the following previous values: yt−7= 28.9, yt−2= 29.2, yt= 29.1. Given these input values, the membership degrees for the fuzzy sets involved in the rule premise are µHIGH1(28.9) = 0.925, µHIGH2(29.2) = 0.990, µHIGH3(29.1) = 0.997, and thus the firing degree of the rule is 0.925, resulting from applying the minimum conjunction operator on the membership degrees of the premise of the rule. This is the firing degree of the rule, and weights the contribution of the rule consequent to the final output of the 20 system. In this case, it is the rule that activates the most, since the three input values clearly fall within the “HIGH” region of the input space. In order to perform an inference, the 15 rules in the rule base are aggregated using the fuzzy mean defuzzification method, i.e., the singleton conclusions of the rules are averaged using their respective firing degrees as weight over the sum of the 15 firing degrees, according to equation 1, yielding ˆyt+1 = 28.728 as prediction of the inference model for this particular data sample of the ENSO series. The significant contribution of the application of an input selection stage based on an effective, nonlinear method can be clearly seen in this case. If we use the DT estimate to select the complexity of the system but no input selection is performed, and thus the inference system has 10 inputs, 3 MFs per input and 98 rules are identified, while the test error is 26% higher and the computational cost increases by more than an order of magnitude. If instead only the three last known values are considered as inputs, 8 MFs per input and 65 rules are identified, while the test error is 31% higher. In the latter case, the use of techniques for finding a better embedding delay [25] and selecting the three inputs accordingly do not provide significant improvements. 4.4 Model Selection Let us now consider the use of the DT estimate for selecting the proper complexity of the tuned fuzzy systems. Rather than using the DT estimate itself, a tolerance band above it is considered. This band is defined by a threshold (DT based threshold, DTBTh) which increases with increasing horizons haccording to equation 2, where DThis the DT estimate for horizon h. DTBTh= (1 + min(0.90,0.15 ·h)) ·DTh(2) 21 For each horizon h, the simplest system that satisfies MSEh≤DTBTh, where MSEh is the training mean square error, is selected as the best autoregressor. This threshold has been defined on the basis of trial and error as a soft limit that favors simplicity to the detriment of accuracy. However, it was found to be robust for all the series analyzed. The definition is based on the following empirical observations: •A tolerance of approximately 15% over the DT estimate for horizon 1 is appropriate. •The best results can be achieved with tolerances increasing with the prediction horizon (particularly for the first 10 predictions approximately). •A tolerance between 80%-100% over the DT estimates provides good results for long-term prediction. We note though that the impact of the threshold is not determining for accuracy (the error increase is of the order of 10-20% at most for any prediction horizon). Similar results can be achieved by selecting a fixed adjustment factor of around 50%-75%. We chose the particular values in equation 2so as to favor model simplicity to the detriment of accuracy. This is further discussed in section 6. For the ENSO series, DTBT1≈1.26 ·10−3and, as shown in figure 5, the fuzzy system with 3 linguistic labels per input is chosen as the best autoregressor for horizon 1. 4.5 Results and Accuracy Comparison Considering now the performance of our methodology for short- and long-term prediction, figure 6shows the normalized DT NNE estimates (NDT-NNE) for prediction horizons up to 50 as well as the training and validation errors of the fuzzy autoregressors built. We note that besides the limitations of the fuzzy modeling techniques being employed, 22 an additional source of error has been introduced in the proposed methodology: the DT based selection of complexity does not guarantee optimal selection under all conditions. Although the fuzzy regressor for horizon 1 prediction that is chosen is the one with the lowest validation error, this is not the case for all horizons. In general, the deviation from the optimal selection depends on the time series being modeled and the prediction horizon. Let us now compare the validation errors of the systems actually selected against the lowest validation errors that could have been achieved for any complexity. This way we can know the order of magnitude of the error due to the imperfection of the DT based complexity selection. Figure 7compares the NDT estimate (a robust estimation of the lowest training error that can be achieved without overfitting), the validation errors of the fuzzy autoregressors selected according to the DT estimate, and the lowest possible validation errors for any number of linguistic labels. Figure 8shows the predictions for the first 50 values after the training set together with a fragment of the actual time series. Finally, we compare the accuracy of fuzzy models against LS-SVM models with the same autoregressor size and input selection. Other models are compared in appendix A. LS-SVMs were built for the same training subset selecting Radial Basis Function (RBF) kernels, grid search as optimization routine and cross-validation as cost function, see [48] for a detailed specification of these and other options. Concrete implementation details will be given in section 6. Figure 9shows the training and generalization errors for LS-SVM and fuzzy models. Averages errors are listed in table 2. Two main conclusions can be drawn from the comparison: •As for generalization capability, the performance of fuzzy autoregressors is clearly better than that of LS-SVM models. There are 4 exceptions: test errors of fuzzy 23 autoregressors are slightly higher (less than 5%) for horizons 14, 18 and 19. However, the overall superiority of fuzzy regressors is clear and specially evident for long-term prediction (beyond horizon 25). •Training and generalization errors are much closer for fuzzy models than for LS-SVM models. For long-term prediction, generalization errors may be even lower than training errors. Also, generalization errors are within approximately 160% of training errors for the worst cases. Thus, training errors of fuzzy models can be trusted as more realistic estimations of the order of the out-of-sample prediction errors. 5 Experimental Results In this section, the proposed concrete implementation of the methodology framework described is applied to a number of varied time series prediction problems from different fields of application, namely the Poland electricity time series prediction benchmark, the monthly averaged sunspot number, the daily averaged aggregated traffic in the Internet2 backbone network, the laser generated data set of the Santa Fe time series competition, and the Mackey-Glass series. In order to ensure reproducibility, the relevant data sources, methods, software tools and parameters used are specified. For every series, models are built to predict the next 50 values. Though one of the major goals of the methodology proposed here is to avoid the need for validation and test series, we will split the series into two subsets in order to assess the out-of-sample prediction performance of the methods being used. Following the methodology illustrated in the previous section, we will summarize the results of the input selection stage and the training and test errors. Results will be compared against those of analogous LS-SVM models built using the same input selection scheme, RBF kernels, grid search as optimization routine and cross-validation as cost function. In appendix A 24 we show a further comparison with other modeling techniques. 5.1 Poland Electricity Benchmark This time series (PolElec henceforward) represents the normalized average daily electricity demand in Poland in the 1990´s. The benchmark consists of a training set of 1400 samples, shown in figure 10(a), and a test set of 201 samples, shown in figure 10(b). It has been shown that the dynamics of this time series is nearly linear [27]. Besides the yearly periodicity, a clear weekly periodicity can be seen on smaller time scales (see figure 10(b)). We will show the results obtained for two different maximum regressor sizes: 7 and 14. In both cases, input selection was performed by exhaustive search of the lowest DT estimate. The number of selected variables is shown in figure 11 Training and test errors of a set of fuzzy autoregressors for horizon 1 are shown in figure 12(a) for different numbers of linguistic labels per input, in the case of a maximum regressor size of 7. The regressor with 5 MFs is selected according to the DT based threshold. Figure 12(b) shows the training and test errors of fuzzy regressors with different numbers of linguistic labels for prediction horizon 7 (also in the case of a maximum regressor size of 7). The system with 2 linguistic labels is selected according to the DT based threshold. However, the system with 3 linguistic labels achieves the lowest test error. This is an illustrative case in which a simpler and less accurate model is selected because of the permissive nature of the DT based threshold. Besides a lower number of linguistic labels, the system with 2 linguistic labels per input has 8 rules, whereas the system with 3 linguistic labels per input has 15 rules. 25 leads us to conclude that in the absence of noise and perturbations, fuzzy inference based autoregression may not be a proper technique if the main objective is approximation accuracy and interpretability is a secondary objective. This type of series is however not common in many real world applications. In addition, the higher accuracy of LS-SVM does not come at no cost. For the MG (9), SFL (10), and SFL (16) series, the construction and optimization of the LS-SVM model requires approximately 37, 35 and 103 times more run time respectively, as shown in table 3. Thus, the methodology proposed here, while clearly less accurate for this kind of series, is still significantly faster and exhibits less overfitting. As far as computational requirements is concerned, the proposed methodology has a very low cost compared against the LS-SVM method. This factor has important practical implications that are often neglected or only partially addressed in the literature. In fact, the high and often unaffordable cost of the LS-SVM and other accurate modeling techniques has recently motivated the development of faster nonlinear learning machines for time series applications [47,20]. A tool, xftsp [38], has been developed that implements the methodology proposed in this paper and provides support for the identification and tuning algorithms included in the Xfuzzy development environment for fuzzy systems [53]. The design of the xftsp tool allows for the use of the wide set of tools available in the Xfuzzy environment for complementary tasks such as visualization, simplification and code generation. This Java based implementation of the methodology presented here is consistently between 1 and 2 orders of magnitude faster than the implementation of LS-SVM used for this study: the optimized C version of the LS-SVMlab1.5 Matlab/C toolbox [30]. Table 3shows the time required to build models with both methods for a subset of the time series considered in this paper. Memory consumption is also much lower for the fuzzy methodology, which enables it to be applied to large training series beyond the few 32 thousand samples current practical limitation of LS-SVM models. The fact that the test results are improved when a DT based threshold higher than the DT estimate itself is introduced, leads us to two remarks on the performance of the identification and tuning stage: •There is likely room for improvement of the identification and tuning procedures. •The DT based threshold can be seen as an aggressiveness index. 1 would be the most aggressive option, most often leading to overfitting and high complexity. Values in the range [1.2,2]DThare reasonable for the identification and learning techniques employed, most often leading to both low complexity and overfitting. The fact that the impact of the DT based threshold is very similar for all the series analyzed leads us to conclude that it is a factor eminently dependent on the identification and learning procedure and its inner limitations. Other methods for fuzzy inference systems identification, tuning and simplification exist, and the ones used in this paper could be improved. This is an area of future research. For complex and noisy time series, it is common that the most simple fuzzy system that can be built (the one with 2 linguistic labels per input) is comparable in accuracy to the LS-SVM model. For example, the fuzzy system with 2 linguistic labels per input for horizon 1 prediction of the PolElec series outperforms LS-SVM with the same input selection. In this case, the test error of the fuzzy regressor is approximately 35% lower. In general, it can be concluded that fuzzy systems with the minimum number of linguistic terms, though not optimal in terms of accuracy, provide a reasonable approximation to the best system that can be built. Thus, it is easy to obtain very simple approximate models that ease the understanding of the time series dynamics. 33 7 Conclusion We have proposed an automatic methodology framework for long-term time series prediction by means of fuzzy inference systems. The use of a nonlinear input selection method yields improved accuracy and interpretability. Experimental results for a concrete implementation of the methodology confirm the satisfactory approximation accuracy and generalization capability of fuzzy regressors. Linguistic interpretability and significantly lower computational requirements are two remarkable advantages over common time series prediction methods. A fundamental advantage of autoregressive time series prediction with fuzzy inference systems is the fact that the models constructed consist of linguistic rules that can be interpreted by humans. For some time series, the most accurate rule bases have a low number of rules (below 10-15 rules), making it easy to draw a linguistic explanation of the system dynamics. Several procedures have been shown to play a key role in achieving good approximation accuracy and low overfitting while keeping the complexity low: variable selection, application of a supervised learning method for tuning after identification, and using DT-NNE for selecting the proper number of linguistic labels per input. Also, when systems have a high number of rules and are thus not interpretable by humans in practice, there is still the possibility to build simpler, approximate models with a degree of accuracy of the same order. The proposed methodology has been shown to clearly outperform LS-SVM based predictions in terms of approximation accuracy except in the particular case of noise-free, stationary and deterministic time series, where fuzzy autoregressors are still significantly faster to build and exhibit less overfitting. 34 A Accuracy Comparison of Different Methods For further comparison, the modeling and prediction accuracy of different modeling techniques for the series analyzed in this paper are shown in table 4. Two alternatives are considered in addition to LS-SVM: OP-ELM and k-NN models. Errors are shown in units relative to that of fuzzy models. The extreme learning machine (ELM) [20] is a simple yet effective learning algorithm for training single-hidden-layer feed-forward artificial neural networks with random hidden nodes. OP-ELM [47] is a methodology based on the ELM, that has been shown to produce models competitive against well-known, accurate techniques, such as LS-SVM and Multilayer Perceptron, while being significantly faster. OP-ELM models were built using the OP-ELM toolbox [33] with the following configuration options: a combination of linear, Gaussian and sigmoid kernels, a maximum of 100 neurons, and data normalization before modeling. k-NN models were generated with 10 as maximum number of neighbors, using the Euclidean distance and 10-fold cross-validation for selecting the best k. LS-SVM models were generated as detailed in sections 5and 6. Considering the out-of-sample or test error, LS-SVM and OP-ELM are in general more accurate than k-NN models, with only a slight exception for the SFL (16) series. LS-SVM are in general more accurate than OP-ELM models, with the clear exceptions of the Sunspots (9) and Sunspots (12) series, where OP-ELM models are the most accurate and also outperform fuzzy models. OP-ELM are also slightly more accurate than LS-SVM for the AbileneI (7) series and slightly less accurate for the AbileneI (14). It should be noted that the Sunspots series is considerably nonstationary and the test set differs from the training set significantly. 35 B Comparison of Different Neuro-Fuzzy Methods In the particular implementation of the proposed methodology used in this paper, the W&M and Levenberg-Marquardt methods have been employed for identification and tuning. This selection has been made on the basis of their satisfactory results in terms of both accuracy and interpretability. Here we compare the performance of these two methods against some other alternatives. Table 5compares the test errors obtained with different methods for identification and tuning of fuzzy inference systems. For easier comparison, errors are shown relative to the errors of the inference models built using the W&M and Levenberg-Marquardt methods within the methodology proposed in this paper. For instance, the 50 fuzzy inference models for the ENSO series, when identified using the W&M method and optimized using the Rprop method, have an average test error 3.8% higher than that shown in section 5, table 2for the W&M and Levenberg-Marquardt methods. Two identification methods are considered in the table: W&M and the method based on subtractive clustering (SC) proposed by Chiu [11]. The W&M method was found to consistently provide better results than other grid partition based methods, such as the algorithm by Higgins and Goodman [18]. The SC method was found to consistently provide better accuracy than other clustering based identification alternatives using the Gath-Geva [1,15], Gustafson-Kessel [17], hard and fuzzy C-means [13] clustering methods. A number of supervised learning algorithms were tested for the tuning substage. For this study, we used the implementations in the Xfuzzy environment, see [39] for a more detailed description of the wide range of methods supported. Among them, we distinguish four classes of methods: gradient descent [32], conjugate gradient, second 36 order or quasi-Newton [3], and algorithms with no derivatives. Table 5shows the test errors for the best option from each of the first three classes of algorithms: Resilient Propagation (Rprop) [42,32], from the gradient descent class, Scaled Conjugated Gradient (SCG) [35], from the conjugate gradient class, and Levenberg-Marquardt (L-M) [3], from the second order class of methods. The following parameters were used for the L-M method: initial Hessian addition 0.1, increase factor 10.0 and decrease factor 0.2. It is worth to mention that none of the statistical or probabilistic methods was found to be competitive in terms of performance, being unable to achieve training errors below the DT based threshold in most cases, within reasonable time bounds. These include the Simulated Annealing method with different cooling schemes, Downhill Simplex and Powell´s methods [39]. We note however that these methods are highly dependent on the values of several parameters that could be explored only partially. From table 5, it is clear that the most accurate models can be obtained with the W&M identification method. Regarding the tuning method, the SCG method is only slightly less accurate than the L-M method, while the Rprop method is slightly more accurate than L-M. The three options have been selected as the most accurate among the set of methods tested, achieving in practice very similar results in terms of accuracy. However, one of the main objectives of this study is to build models that are accurate, yet as simple as possible. If we look at the number of MFs required to achieve such degrees of accuracy, the L-M method is more efficient, as can be seen in table 6. The table shows three measures of complexity of the fuzzy systems generated with each method: the number of MFs, as well as the percent shift of the centers and the percent change of widths of the MFs after tuning, with respect to the initial uniform grid partition. For a lower number of input MFs, higher changes in the shapes of the MFs can be expected. As a conclusion, we used the L-M method in this work because it produces systems with 37 the lowest number of MFs, while being almost as accurate as possible with the methods tested. We note however that the Rprop method is almost equivalent in terms of the number of MFs required. Rprop would be thus a good alternative, yielding models that are slightly more complex and accurate. Acknowledgements The first author is supported by a Marie Curie Intra-European Fellowship for Career Development (grant agreement PIEF-GA-2009-237450) within the European Community´s Seventh Framework Programme (FP7/20072013). Most of this work was done while the first author was with the Microelectronics Institute of Seville, IMSE-CNM, CSIC – Scientific Research Council. C. Americo Vespucio s/n. Parque Tecnologico Cartuja. E-41092 Seville, Spain. This work has been supported in part by project TEC2008-04920/MICINN from the Spanish Ministry of Science and Innovation, as well as project P08-TIC-03674 and grants IAC07-I-0205:33080 and IAC08-II-3347:56263 from the Andalusian regional Government. 38 References [1] J. Abonyi, R. Babuska, F. Szeifert, Modified Gath-Geva Fuzzy Clustering for Identification of Takagi-Sugeno Fuzzy Models, IEEE Transactions on Systems, Man and Cybernetics, Part B 32 (5) (2002) 612–621. [2] P. P. Angelov, D. P. Filev, An approach to online identification of takagi-sugeno fuzzy models, IEEE Transactions on Systems, Man and Cybernetics - Part B: Cybernetics 34 (1) (2004) 484–498. [3] R. Battiti, First and Second Order Methods for Learning: Between Steepest Descent and Newton’s Method, Neural Computation 4 (2) (1992) 141–166. [4] I. Baturone, F. J. Moreno-Velo, A. Gersnoviez, A CAD Approach to Simplify Fuzzy System Descriptions, in: 15th IEEE International Conference on Fuzzy Systems (FUZZ-IEEE’06), 2006. [5] R. E. Bellman, Dynamic Programming, 1957th ed., Republished, Dover Publications Inc., 2003, ISBN: 0486428095. [6] G. Box, G. M. Jenkins, G. Reinsel, Time Series Analysis: Forecasting & Control, Prentice Hall; 3rd edition, 1994, ISBN: 0130607746. [7] J. Casillas, O. Cord´ on, F. Herrera, L. Magdalena (eds.), Interpretability Issues in Fuzzy Modeling, Studies in Fuzziness and Soft Computing, Springer Verlag, Berlin, Germany, 2003, ISBN: 978-3-540-02932-8. [8] Y.-H. O. Chang, B. M. Ayyub, Fuzzy regression methods - a comparative assessment, Fuzzy Sets and Systems 119 (2) (2001) 187–203. [9] C. Chatfield, The Analysis of Time Series. An Introduction, CRC Press, 2003, Sixth edition, ISBN: 1-58488-317-0. [10] C. Chen, L.-M. Liu, Joint Estimation of Model Parameters and Outlier Effects in Time Series, Journal of the American Statistical Association 88 (421) (1993) 284–297. 39 [11] S. L. Chiu, A Cluster Estimation Method with Extension to Fuzzy Model Identification, in: IEEE Conference on Fuzzy Systems, 1994. IEEE World Congress on Computational Intelligence, Orlando, FL, USA, 1994. [12] S. L. Chiu, Fuzzy Model Identification Based on Cluster Estimation, Journal of Intelligent & Fuzzy Systems 2 (3) (1994) 267–278. [13] J. V. de Oliveira, W. Pedrycz (eds.), Advances in Fuzzy Clustering and its Applications, John Wiley & Sons, Ltd., West Sussex, England, 2007, ISBN: 978-0-470-02760-8. [14] ESTSP´07 European Symposium on Time Series Prediction: Prediction Competition (Mar. 2008). URL http://www.estsp.org [15] I. Gath, A. B. Geva, Unsupervised Optimal Fuzzy Clustering, IEEE Transactions on Pattern Analysis and Machine Intelligence 11 (7) (1989) 773–780. [16] S. Guillaume, Designing fuzzy inference systems from data: An interpretability-oriented review, IEEE Transactions on Fuzzy Systems 9 (3) (2001) 426–443. [17] E. E. Gustafson, W. C. Kessel, Fuzzy Clustering with a Fuzzy Covariance Matrix, in: 17th Symposium on Adaptive Processes, 1978 IEEE Conference on Decision and Control , San Diego, CA, 1978. [18] C. M. Higgins, R. M. Goodman, Fuzzy Rule-Based Networks for Control, IEEE Transactions on Fuzzy Systems 2 (1) (1994) 82–88. [19] T.-P. Hong, C.-Y. Lee, Induction of fuzzy rules and membership functions from training examples, Fuzzy Sets and Systems 84 (1) (1996) 33–47. [20] G.-B. Huang, Q. Y. Zhu, C. K. Siew, Extreme learning machine: Theory and applications, Neurocomputing 70 (1–3) (2006) 489–501. [21] The Internet2 Observatory (Jul. 2008). URL http://www.internet2.edu/observatory/ 40 [22] J.-S. R. Jang, ANFIS: Adaptive-Network-Based Fuzzy Inference System, IEEE Transactions on Systems, Man and Cybernetics 23 (3) (1993) 665–685. [23] J.-S. R. Jang, C.-T. Sun, E. Mizutani, Neuro-Fuzzy and Soft Computing A Computational Approach to Learning and Machine Intelligence, Prentice Hall, Upper Saddle River, New Jersey, 1997, ISBN 0-13-261066-3. [24] A. J. Jones, New Tools in Non-linear Modelling and Prediction, Computational Management Science 2 (1) (2004) 109–149. [25] H. Kantz, T. Schreiber, Nonlinear Time Series Analysis, 2nd ed., Cambridge University Press, Cambridge, UK, 2004. [26] N. K. Kasabov, Q. Song, DENFIS: Dynamic Evolving Neural-Fuzzy Inference System and Its Application for Time-Series Prediction, IEEE Transactions on Fuzzy Systems 10 (2) (2002) 144–154. [27] A. Lendasse, J. Lee, V. Wertz, M. Verleyssen, Forecasting Electricity Consumption using Nonlinear Projection and Self-Organizing Maps, Neurocomputing 48 (1) (2002) 299–311. [28] E. Liiti¨ ainen, A. Lendasse, F. Corona, Non-parametric Residual Variance Estimation in Supervised Learning, in: IWANN 2007, International Work-Conference on Artificial Neural Networks, San Sebasti´ an, Spain, 2007. [29] E. Liiti¨ ainen, A. Lendasse, F. Corona, Bounds on the mean power-weighted nearest neighbour distance, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 464 (2097) (2008) 2293–2301. [30] Least Squares - Support Vector Machines Matlab/C Toolbox (Apr. 2008). URL http://www.esat.kuleuven.ac.be/sista/lssvmlab [31] M. C. Mackey, L. Glass, Oscillations and Chaos in Physiological Control Systems, Science 197 (4300) (1977) 287–289. 41 Table 3 Run time (in seconds) required to build models for prediction horizons 1-50. All tests were run on the same system, with no significant competing load. Maximum regressor size specified between parenthesis. Series LS-SVMlab1.5 Fuzzy inference ENSO (10) 3.45 ·1051.05 ·104 PolElec (7) 3.04 ·1051.05 ·104 PolElec (14) 9.91 ·1052.30 ·104 Sunspots (9) 3.10 ·1051.04 ·104 Sunspots (12) 2.42 ·1051.22 ·104 AbileneI (7) 1.40 ·1051.75 ·103 AbileneI (12) 1.27 ·1054.69 ·103 SFL (10) 1.28 ·1063.49 ·104 SFL (16) 1.61 ·1064.55 ·104 MG (9) 3.64 ·1053.54 ·103 48 Table 4 Accuracy comparison of different nonlinear modeling techniques. Training and test errors, averaged for horizons 1 through 50, of different nonlinear models. Training and test errors are expressed relative to the training and test errors, respectively, of fuzzy models built using the W&M and Levenberg-Marquardt methods within the proposed methodology. Absolute errors for the fuzzy and LS-SVM models were given in section 5, table 2. LS-SVM OP-ELM k-NN Series Training Test Training Test Training Test ENSO (10) 0.41 1.56 0.76 1.91 0.18 2.30 PolElec (7) 0.68 1.99 0.89 2.52 0.39 2.87 PolElec (14) 0.66 1.78 0.71 2.29 0.34 2.86 Sunspots (9) 0.79 1.25 0.92 0.88 0.22 1.34 Sunspots (12) 0.61 1.19 0.90 0.93 0.13 1.66 AbileneI (7) 0.59 1.43 0.74 1.39 0.44 2.32 AbileneI (12) 0.55 1.43 0.85 1.51 0.29 2.06 SFL (10) 0.15 0.51 1.23 1.14 0.27 1.26 SFL (16) 0.06 0.45 1.56 1.18 0.22 1.16 MG (9) 0.06 0.21 0.25 1.52 0.19 1.63 49 Table 5 Accuracy comparison of different methods for building fuzzy inference models. Test errors, averaged for horizons 1 through 50, are expressed relative to the test errors of the equivalent fuzzy models built using the W&M and Levenberg-Marquardt methods within the proposed methodology. Results are shown for two identification methods: W&M and Subtractive Clustering (SC) based, as well as three tuning methods: Levenberg-Marquardt (L-M), Scaled Conjugated Gradient (SCG), and Resilient Propagation (Rprop). W&M SC SCG Rprop L-M SCG Rprop ENSO (10) 1.003 1.038 1.120 1.065 1.249 PolElec (7) 0.958 0.991 1.035 1.062 1.058 PolElec (14) 0.941 0.964 1.045 1.085 1.070 Sunspots (9) 1.038 0.997 1.139 1.132 1.148 Sunspots (12) 0.996 0.953 1.068 1.100 1.061 AbileneI (7) 1.004 0.922 2.258 2.739 1.080 AbileneI (12) 1.157 1.112 3.177 3.092 1.381 SFL (10) 1.000 0.988 1.055 1.169 1.102 SFL (16) 0.995 0.982 1.089 1.183 1.215 MG (9) 1.051 1.011 1.082 1.184 1.123 Average 1.014 0.997 1.407 1.481 1.153 50 Table 6 Complexity comparison of different methods for building fuzzy inference models. Three supervised learning methods are compared: Levenberg-Marquardt (L-M), Scaled Conjugated Gradient (SCG), and Resilient Propagation (Rprop). Three measures of complexity and interpretability are shown: number of membership functions (#MFs), percent center shift (∆cj,k), and percent width change (∆σj,k). L-M SCG Rprop #MF ∆cj,k ∆σj,k #MF ∆cj,k ∆σj,k #MF ∆cj,k ∆σj,k ENSO (10) 3.84 8.53 21.7 4.26 6.35 21.2 4.18 7.40 18.0 PolElec (7) 3.08 12.3 13.2 3.14 6.11 8.82 3.10 11.0 10.3 PolElec (14) 3.10 12.8 15.8 3.08 5.38 9.12 3.20 11.4 10.7 Sunspots (9) 3.02 7.40 9.72 3.02 4.66 7.14 3.00 5.89 9.72 Sunspots (12) 3.03 9.20 12.1 3.03 5.64 8.04 3.02 6.25 10.2 AbileneI (7) 3.13 5.87 9.72 3.40 4.02 8.88 3.12 6.60 8.94 AbileneI (12) 4.60 5.13 10.9 5.26 4.99 15.4 4.44 6.73 15.2 SFL (10) 5.03 6.03 26.5 5.94 5.26 27.8 4.97 6.81 22.7 SFL (16) 4.77 5.88 26.1 5.36 5.31 26.4 4.71 6.57 22.4 MG (9) 3.72 12.1 17.7 4.34 7.40 16.0 4.06 10.8 16.0 Average 3.732 8.52 16.3 4.832 5.53 14.89 3.784 7.95 14.42 51 List of Figures 1 Methodology Framework for Time Series Prediction. 54 2 ESTSP´07 competition data set (ENSO series, 875 samples). 55 3 ENSO: Number of selected variables for horizon up to 50. DT based selection with exhaustive search. Maximum regressor size 10. 56 4 ENSO: Membership functions for the three inputs and the output of the autoregressor for prediction horizon 1, after the tuning substage. The resulting input functions are approximately the same as the initial uniform partition defined before the tuning substage. 57 5 ENSO: Errors for horizon 1, exhaustive DT based selection of inputs. Continuous line: training error. Dashed line: validation error. 58 6 ENSO: NDT estimates (∗), training (+) and validation (×) errors of fuzzy autoregressors. Maximum regressor size 10. DT based selection of inputs. 59 7 ENSO: NDT estimates (∗), test errors for the selected fuzzy autoregressors (+), validation errors for the optimal complexity selections (x). Maximum regressor size 10. DT based selection of inputs. 60 8 ENSO: Prediction of 50 values after the training set. Continuous line (+): actual time series. Dashed line (×): predictions. 61 9 ENSO: comparison of the proposed methodology against LS-SVM. Generalization errors of LS-SVM models (+). Generalization errors of fuzzy models (). Training errors of fuzzy models (∗). Training errors of LS-SVM models (×). 62 10 PolElec: training and test series. 63 11 PolElec: number of selected variables (exhaustive DT based selection). Continuous line: regressor size 7. Dashed line: regressor size 14. 64 12 PolElec: training (continuous line) and test (dashed line) errors against linguistic labels per input. Exhaustive DT based selection of variables with maximum regressor size 7. 65 13 PolElec: NDT estimates (∗), training (+) and test (×) errors of fuzzy autoregressors. Exhaustive DT based selection of inputs. 66 14 Sunspots: training (first 1000 samples) and test (last 2098 samples) series. 67 52 15 Sunspots: number of selected variables (exhaustive DT based selection). Continuous line: regressor size 9. Dashed line: regressor size 12. 68 16 Sunspots: NDT estimates (∗), training (+) and test (×) errors of fuzzy autoregressors. Exhaustive DT based selection of inputs. 69 17 Internet2: daily averaged aggregated incoming traffic in the Abilene backbone for 1458 days. Training series (first 728 values) and test series (last 730 values). 70 18 Internet2: number of selected variables (exhaustive DT based selection). Continuous line: regressor size 7. Dashed line: regressor size 12. 71 19 Internet2: NDT estimates (∗), training (+) and test (×) errors of fuzzy autoregressors. Exhaustive DT based selection of inputs. 72 20 SFL: training and test series. 73 21 SFL: Number of selected variables. Continuous line: exhaustive DT search with maximum regressor size 10. Dashed line: forward-backward DT search with maximum regressor size 16. 74 22 SFL: NDT estimates (∗), training (+) and test (×) errors of fuzzy autoregressors. Exhaustive DT based selection of inputs for a maximum regressor size of 10. For regressor size 16, the exhaustive search for size 10 is extended using a forward-backward DT based search. 75 23 MG: fragment of the Mackey-Glass series (1500 samples). The first 500 samples are selected as training set. The remaining 1000 samples are selected as test set. 76 24 MG: Number of selected variables for horizons up to 50. Exhaustive DT based selection of inputs. Maximum regressor size 9. 77 25 MG: NDT estimates (∗), training (+) and test (×) errors of fuzzy autoregressors. Exhaustive DT based selection of inputs. Maximum regressor size 9. 78 53 Fig. 1. Methodology Framework for Time Series Prediction. 54 18 20 22 24 26 28 30 0 100 200 300 400 500 600 700 800 900 Temperature Week Fig. 2. ESTSP´07 competition data set (ENSO series, 875 samples). 55 1 2 3 4 5 6 7 8 9 10 5 10 15 20 25 30 35 40 45 50 Selected variables Prediction horizon Fig. 3. ENSO: Number of selected variables for horizon up to 50. DT based selection with exhaustive search. Maximum regressor size 10. 56 (a) yt−7(b) yt−2(c) yt(d) yt+1 Fig. 4. ENSO: Membership functions for the three inputs and the output of the autoregressor for prediction horizon 1, after the tuning substage. The resulting input functions are approximately the same as the initial uniform partition defined before the tuning substage. 57 1 2 3 4 5 6 7 8 9 10 11 12 13 14 5 10 15 20 25 30 35 40 45 50 Selected variables Prediction horizon Fig. 11. PolElec: number of selected variables (exhaustive DT based selection). Continuous line: regressor size 7. Dashed line: regressor size 14. 64 0 0.001 0.002 0.003 0.004 0.005 0.006 2 3 4 5 6 7 8 9 10 11 12 13 14 15 NMSE Linguistic labels per input variable (a) Horizon 1 0.004 0.005 0.006 0.007 0.008 2 3 4 5 6 7 8 9 10 11 NMSE Linguistic labels per input variable (b) Horizon 7 Fig. 12. PolElec: training (continuous line) and test (dashed line) errors against linguistic labels per input. Exhaustive DT based selection of variables with maximum regressor size 7. 65 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 5 10 15 20 25 30 35 40 45 50 NMSE Prediction horizon (a) Maximum regressor size 7 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 5 10 15 20 25 30 35 40 45 50 NMSE Prediction horizon (b) Maximum regressor size 14 Fig. 13. PolElec: NDT estimates (∗), training (+) and test (×) errors of fuzzy autoregressors. Exhaustive DT based selection of inputs. 66 0 50 100 150 200 250 0 100 200 300 400 500 600 700 800 900 1000 Sunspot number Months Fig. 14. Sunspots: training (first 1000 samples) and test (last 2098 samples) series. 67 1 2 3 4 5 6 7 8 9 10 11 12 5 10 15 20 25 30 35 40 45 50 Selected variables Prediction horizon Fig. 15. Sunspots: number of selected variables (exhaustive DT based selection). Continuous line: regressor size 9. Dashed line: regressor size 12. 68 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 5 10 15 20 25 30 35 40 45 50 NMSE Prediction horizon (a) Maximum regressor size 9 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 5 10 15 20 25 30 35 40 45 50 NMSE Prediction horizon (b) Maximum regressor size 12 Fig. 16. Sunspots: NDT estimates (∗), training (+) and test (×) errors of fuzzy autoregressors. Exhaustive DT based selection of inputs. 69 0 2e+08 4e+08 6e+08 8e+08 1e+09 1.2e+09 1.4e+09 1.6e+09 1.8e+09 2e+09 2.2e+09 0 250 500 750 1000 1250 Aggregated incoming traffic (bps) Day Fig. 17. Internet2: daily averaged aggregated incoming traffic in the Abilene backbone for 1458 days. Training series (first 728 values) and test series (last 730 values). 70 1 2 3 4 5 6 7 8 9 10 11 12 5 10 15 20 25 30 35 40 45 50 Selected variables Prediction horizon Fig. 18. Internet2: number of selected variables (exhaustive DT based selection). Continuous line: regressor size 7. Dashed line: regressor size 12. 71 0 0.005 0.01 0.015 0.02 0.025 0.03 5 10 15 20 25 30 35 40 45 50 NMSE Prediction horizon (a) Maximum regressor size 7 0 0.005 0.01 0.015 0.02 0.025 5 10 15 20 25 30 35 40 45 50 NMSE Prediction horizon (b) Maximum regressor size 12 Fig. 19. Internet2: NDT estimates (∗), training (+) and test (×) errors of fuzzy autoregressors. Exhaustive DT based selection of inputs. 72 0 50 100 150 200 250 300 0 100 200 300 400 500 600 700 800 900 1000 Intensity time (a) Training series (1000 samples) 0 50 100 150 200 250 300 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 Intensity time (b) Test series (9093 samples) Fig. 20. SFL: training and test series. 73