An ARMA Type Fuzzy Time Series Forecasting Method Based on Particle Swarm Optimization

Abstract

In the literature, fuzzy time series forecasting models generally include fuzzy lagged variables. Thus, these fuzzy time series models have only autoregressive structure. Using such fuzzy time series models can cause modeling error and bad forecasting performance like in conventional time series analysis. To overcome these problems, a new first-order fuzzy time series which forecasting approach including both autoregressive and moving average structures is proposed in this study. Also, the proposed model is a time invariant model and based on particle swarm optimization heuristic. To show the applicability of the proposed approach, some methods were applied to five time series which were also forecasted using the proposed method. Then, the obtained results were compared to those obtained from other methods available in the literature. It was observed that the most accurate forecast was obtained when the proposed approach was employed.

1. Introduction

Fuzzy time series firstly proposed by Song and Chissom [1] can be divided into two subclasses time variant and time invariant. Fuzzy time series method generally embodies three stages such as fuzzification, determination of fuzzy relations, and defuzzification stages. In the fuzzification stage, observations of time series are fuzzified. Fuzzy relations between the observations are defined in the stage of determination of fuzzy relations. Finally, the calculated fuzzy forecasts are defuzzified in the defuzzification phase. If one or more of these stages can be improved, the performance of the method will increase. Therefore, new fuzzy time series approaches have been proposed by making contributions to these stages in the literature.

In the literature, various methods have been proposed to fuzzify observations. While fixed interval lengths are used in Song and Chissom [1–3], Chen [4], Huarng [5], Chen [6], Tsaur et al. [7], Singh [8], and Egrioglu et al. [9, 10], dynamic length of interval lengths is employed in Huarng and Yu [11], Davari et al. [12], Yolcu et al. [13], Kuo et al. [14, 15], Park et al. [16], Hsu et al. [17], and Huang et al. [18] in order to partition the universe of discourse. Also, Cheng et al. [19], Li et al. [20], Egrioglu et al. [21], Chen and Tanuwijaya [22], and Bang and Lee [23] used some methods based on clustering algorithms.

Song and Chissom [3] exploited feed forward neural networks to defuzzify fuzzy forecasts. Chen [4], Huarng [5], and Huarng and Yu [11] utilized centroid method in the defuzzification stage. Besides, Cheng et al. [24] and Aladag et al. [25] used a different technique based on adaptive expectation and centroid methods for fuzzification.

Establishing fuzzy relations plays important role in the forecasting performance of the method. In this phase, Song and Chissom [1] utilized fuzzy relationship matrix, and Sullivan and Woodall [26] used transition matrices based on Markov chain instead of fuzzy relationship matrix. Chen [4] suggested an approach in which fuzzy logic group relationship tables are employed to define fuzzy relations, and Cheng et al. [24], Huarng [5], Huarng and Yu [11], Yu [27], and Egrioglu et al. [21] also used fuzzy logic group relationship tables. Huarng and Yu [28] and Aladag et al. [29] preferred to utilize feed forward artificial neural networks in this stage. Aladag et al. [30] used a different type of artificial neural networks which is Elman recurrent neural networks. In the determination of fuzzy relations, Yu and Huarng [31] and Yolcu et al. [32] proposed different approaches in which feed forward artificial neural networks using membership values are used. Eğrioğlu [33] utilized genetic algorithms to define fuzzy relations. Moreover, many soft computing techniques have been used for forecasting in the literature. Yang et al. [34], and Yang et al. [35], and Yang et al. [36] are some of them.

In all of the studies mentioned previously, fuzzy time series forecasting models have autoregressive (AR) structure. However, many real life fuzzy time series can contain both autoregressive and moving average (MA) structures. Using forecasting models that have only AR structure for these time series can produce insufficient results. To analyze such time series, it is needed to use an ARMA type fuzzy time series forecasting model that includes both AR and MA structures. In this study, a novel fuzzy time series forecasting approach based on ARMA type fuzzy time series forecasting model is proposed to increase forecasting accuracy. Some fuzzy time series approaches in which fuzzy logic group relationship tables are employed disregard the fuzzy set theory since fuzzy sets’ elements whose membership value is 1 are only taken into account in the fuzzification phase. Therefore, information loss occurs, and it is expected to decrease the explanation power of the model. In order to overcome these problems, the particle swarm optimization method is employed in the proposed approach to establish fuzzy relations. In addition, the computational cost of the proposed approach is very reasonable since it does not need to perform complex matrix operations when the particle swarm optimization method is used. Another advantage of the proposed method is that it uses the fuzzy 𝑐-means clustering method in the fuzzification phase. In the fuzzification phase of the fuzzy time series approach, there are some common problems such as the decision of the number of intervals, arbitrary determination of interval length, and arbitrary choice of degrees of membership. Since the fuzzy 𝑐-means clustering technique is employed in the proposed approach, it works in a systematic way.

In the real life, some time series have long-range dependence. Long-range-dependent time series is forecasted by using special linear models. ARFIMA models, which are used to forecast long-range-dependent time series, were firstly introduced by Granger and Joyeux [37]. General properties of ARFIMA models were given by Hosking [38, 39] and Beran [40]. First studies are concerned with estimation of fractional differencing parameter in fractional white noise processes. 𝑅/𝑆 statistic was proposed in Hurst [41]. Egrioglu and Gunay [42], and Li [43], Li and Zhao [44] are important papers about fractional processes. In this study, time series are examined by using 𝑅/𝑆 hypothesis test to determine long-range dependence in the application. The performance of proposed method is examined according to long-range dependence.

The next section presents basic definitions of fuzzy time series and fuzzy-𝑐 means clustering method. In Section 3, the modified particle swarm optimization algorithm is given. Section 4 introduces the proposed approach based on ARMA type fuzzy time series forecasting model and particle swarm optimization. The implementation is given in Section 5. Finally, the last section concludes the paper.

2. Fuzzy Time Series Basic Definitions and Fuzzy 𝐶-Means Clustering

The fuzzy time series was firstly introduced by Song and Chissom [1]. The fundamental definitions of fuzzy time series, time variant, and time invariant fuzzy time series definitions are presented in what follows.

Definition 1. Let 𝑌(𝑡)(𝑡=…,0,1,2,…), a subset of real numbers, be the universe of discourse on which fuzzy sets 𝑓_𝑗(𝑡) are defined. If 𝐹(𝑡) is a collection of 𝑓₁(𝑡),𝑓₂(𝑡),…, then 𝐹(𝑡) is called a fuzzy time series defined on 𝑌(𝑡).

Definition 2. Suppose that 𝐹(𝑡) is caused by 𝐹(𝑡−1) only; that is, 𝐹(𝑡−1)→𝐹(𝑡). Then, this relation can be expressed as 𝐹(𝑡)=𝐹(𝑡−1)∘𝑅(𝑡,𝑡−1), where 𝑅(𝑡,𝑡−1) is the fuzzy relationship between 𝐹(𝑡−1) and 𝐹(𝑡) and 𝐹(𝑡)=𝐹(𝑡−1)∘𝑅(𝑡,𝑡−1) is called the first-order model of 𝐹(𝑡). “∘” represents max-min composition of fuzzy sets.

Definition 3. Suppose that 𝑅(𝑡,𝑡−1) is a first-order model of 𝐹(𝑡). If for any 𝑡,𝑅(𝑡,𝑡−1) is independent of ; that is, for any 𝑡,𝑅(𝑡,𝑡−1)=𝑅(𝑡−1,𝑡−2); then 𝐹(𝑡) is called a time invariant fuzzy time series; otherwise, it is called a time variant fuzzy time series.

Song and Chissom [2] firstly introduced an algorithm based on the first-order model for forecasting time invariant 𝐹(𝑡). In Song and Chissom [2], the fuzzy relationship matrix 𝑅(𝑡,𝑡−1)=𝑅 is obtained by many matrix operations. The fuzzy forecasts are obtained based on max-min composition as in what follows: 𝐹(𝑡)=𝐹(𝑡−1)∘𝑅. (1)

The dimension of 𝑅 matrix depends on the number of fuzzy sets. The number of fuzzy sets equals the number of intervals that is composed of universe of discourse. The more fuzzy sets are used, the more matrix operations are needed to obtain 𝑅 matrix. When the number of fuzzy set is high, using the method proposed by Song and Chissom [2] considerably increases the computational cost.

For fuzzification, partition of universe of discourse method is used in the method proposed by Song and Chissom [1]. However, there are several problems related to the decomposition of universe of discourse. These problems are the determination of the number of intervals, arbitrarily choice of interval length, and membership degrees. In order to deal with these problems, Cheng et al. [19] and Li et al. [20] used fuzzy 𝑐-means clustering method for fuzzification. The fuzzy 𝑐-means clustering method was firstly introduced by Bezdek [45]. This clustering method is the most widely used one. In this method, fuzzy clustering is conducted by minimizing the least squared errors within groups. Let 𝑢_𝑖𝑗 be the membership values, 𝑣_𝑖 the center of cluster, 𝑛 the number of variables, and 𝑐 the number of clusters. Then, the objective function, which is minimized in fuzzy clustering, is 𝐽_𝛽(𝑋,𝑉,𝑈)=^𝑐𝑛_𝑖=1𝑗=1𝑢^𝛽_𝑖𝑗𝑑²𝑥_𝑗,𝑣_𝑖, (2) where 𝛽 is a constant (𝛽>1) and called the fuzzy index. 𝑑(𝑥_𝑗,𝑣_𝑖) is a similarity measure between an observation and the center of corresponding fuzzy cluster. The objective function 𝐽_𝛽 is minimized subject to constraints given in what follows: 0≤𝑢_𝑖𝑗≤1,∀𝑖,𝑗,0<^𝑛_𝑗=1𝑢_𝑖𝑗≤𝑛,∀𝑖,_𝑐𝑖=1𝑢_𝑖𝑗=1,∀𝑗. (3)

In fuzzy 𝑐-means clustering method, to solve the minimization problem given previously, an iterative algorithm is used. In each iteration, the values of 𝑣_𝑖 and 𝑢_𝑖𝑗 are updated by using the formulas given in (4) and (5), respectively: 𝑣_𝑖=∑^𝑛_𝑗=1𝑢^𝛽_𝑖𝑗𝑥_𝑗∑^𝑛_𝑗=1𝑢^𝛽_𝑖𝑗𝑢_𝑖𝑗,=1∑^𝑐_𝑘=1𝑑𝑥_𝑗,𝑣_𝑖/𝑑𝑥_𝑗,𝑣_𝑘^2/(𝛽−1). (4)(5)

3. The Particle Swarm Optimization

Particle swarm optimization, which is a population-based heuristic algorithm, was firstly proposed by Kennedy and Eberhart [46]. Distinguishing feature of this heuristic algorithm is that it simultaneously examines different points in different regions of the solution space to find the global optimum solution. Local optimum traps can be avoided because of this feature.

In a few fuzzy time series studies, particle swarm optimization method has been exploited in fuzzification phase. While the particle swarm optimization method was employed by Davari et al. [12] for fuzzification in the first-order fuzzy time series forecasting model, Kuo et al. [14] utilized the method in high-order models. In the fuzzification phase, Kuo et al. [15] utilized the particle swarm optimization method in both the first- and the high-order models. Park et al. [16] used the same method for fuzzification in a two-factor high-order model. However, the particle swarm optimization method has never been used to establish fuzzy logic relations in the literature. In time variant fuzzy time series forecasting method proposed in this study, the modified particle swarm optimization whose algorithm is given later is exploited. The modified particle swarm optimization algorithm has time varying inertia weight like in Shi and Eberhart [47]. In a similar way, this algorithm also has time varying acceleration coefficient like in Ma et al. [48].

Algorithm 4. The algorithm of the modified particle swarm optimization.

Step 1. Positions of each 𝑘th (𝑘=1,2,…,𝑝𝑛) particles’ positions are randomly determined and kept in a vector 𝑋_𝑘 given as follows: 𝑋_𝑘=𝑥_𝑘,1,𝑥_𝑘,2,…,𝑥_𝑘,𝑑,𝑘=1,2,…,𝑝𝑛, (6) where 𝑥^𝑘_𝑖(𝑖=1,2,…,𝑑) represents 𝑖th position of 𝑘th particle. 𝑝𝑛 and 𝑑 represent the number of particles in a swarm and positions, respectively.

Step 2. Velocities are randomly determined and stored in a vector 𝑉_𝑘 given in what follows: 𝑉_𝑘=𝑣_𝑘,1,𝑣_𝑘,2,…,𝑣_𝑘,𝑑,𝑘=1,2,…,𝑝𝑛. (7)

Step 3. According to the evaluation function, 𝑃best and 𝐺best particles given in (8) and (9), respectively, are determined: 𝑃best_𝑘=𝑝_𝑘,1,𝑝_𝑘,2,…,𝑝_𝑘,𝑑𝑝,𝑘=1,2,…,𝑝𝑛,𝐺best=_𝑔,1,𝑝_𝑔,2,…,𝑝_𝑔,𝑑, (8)(9) where 𝑃best_𝑘 is a vector that stores the positions corresponding to the 𝑘th particle’s best individual performance and 𝐺best represents the best particle, which has the best evaluation function value, found so far.

Step 4. Let 𝑐₁ and 𝑐₂ represent cognitive and social coefficients, respectively, and 𝑤 is the inertia parameter. Let (𝑐_1𝑖,𝑐_1𝑓),(𝑐_2𝑖,𝑐_2𝑓), and (𝑤₁,𝑤₂) be the intervals which include possible values for 𝑐₁, 𝑐₂, and 𝑤, respectively. At each iteration, these parameters are calculated by using the formulas given in (10): 𝑐₁=𝑐_1𝑓−𝑐_1𝑖𝑡max𝑡+𝑐_1𝑖,𝑐₂=𝑐_2𝑓−𝑐_2𝑖𝑡max𝑡+𝑐_2𝑖,𝑤=𝑤₂−𝑤₁max𝑡−𝑡max𝑡+𝑤₁, (10) where max𝑡 and 𝑡 represent maximum iteration number and current iteration number, respectively.

Step 5. Values of velocities and positions are updated by using the formulas given in (11) and (12), respectively. 𝑣^𝑡+1_𝑖,𝑑=𝑤×𝑣^𝑡_𝑖,𝑑+𝑐₁×rand₁×𝑝_𝑖,𝑑−𝑥_𝑖,𝑑+𝑐₂×rand₂×𝑝_𝑔,𝑑−𝑥_𝑖,𝑑,𝑥^𝑡+1_𝑖,𝑑=𝑥_𝑖,𝑑+𝑣^𝑡+1_𝑖,𝑑, (11)(12) where rand₁and rand₂ are random values from the interval [01].

Step 6. Steps 3 to 5 are repeated until a predetermined maximum iteration number (max𝑡) is reached.

We would like to note that the aim of this study is not to propose a new particle swarm optimization algorithm. In the literature, it was shown that using some time varying parameters can increase the convergence speed of the algorithm (Shi and Eberhart [47]; Ma et al. [48]). Therefore, we just used both the time varying acceleration coefficient (Ma et al. [48]) and the time varying inertia weight (Shi and Eberhart [47]) in the standard particle swarm optimization as mentioned previously. Then, this modified particle swarm optimization method was utilized to calculate membership values in the fuzzy relationship matrix as addressed in the next section. It is not claimed in the paper that a better particle swarm optimization algorithm was improved.

4. The Proposed Approach

The fuzzy time series forecasting methods in the literature have only AR structure that includes lagged fuzzy variables. On the other hand, it will be wiser to use a model which includes both AR and MA structures in order to forecast real life time series that contain both structures. In the literature, we firstly define a new model which includes both AR and MA terms in Definition 5.

Definition 5. Let fuzzy time series 𝐹(𝑡) be caused by lagged fuzzy time series 𝐹(𝑡−1) and lagged absolute fuzzy error series 𝜀(𝑡−1); then the relationship can be expressed as 𝐹(𝑡−1),𝜀(𝑡−1)⟶𝐹(𝑡). (13) This model given in (13) is called an ARMA(1,1) fuzzy time series forecasting model. Thus, definition of fuzzy relations for this model is presented in Definition 6.

Definition 6. Suppose that 𝐹(𝑡) is caused by only 𝐹(𝑡−1) and 𝜀(𝑡−1); that is, 𝐹(𝑡−1), 𝜀(𝑡−1)→𝐹(𝑡). Then, this relation can be expressed as 𝐹(𝑡)=𝐹(𝑡−1)∘𝑅₁(𝑡,𝑡−1)∩𝜀(𝑡−1)∘𝑅₂(𝑡,𝑡−1), (14) where 𝑅₁(𝑡,𝑡−1)(𝑅₁=[𝑟_𝑖𝑗1],𝑖=1,2,…,𝑐,𝑗=1,2,…,𝑐) is the fuzzy relationship matrix between 𝐹(𝑡−1) and 𝐹(𝑡) and 𝑅₂(𝑡,𝑡−1)(𝑅₂=[𝑟_𝑖𝑗2],𝑖=1,2,…,𝑐,𝑗=1,2,…,𝑐) is the fuzzy relationship matrix between 𝜀(𝑡−1) and 𝐹(𝑡).

In this study, a new fuzzy time series forecasting approach is suggested to forecast ARMA(1,1) fuzzy time series forecasting model in (13). The proposed approach provides some advantages which can be given as follows.(i)Since the model given in (13) is employed in the proposed approach, it is more proper to use this method for forecasting real life time series that contain both AR and MA structures.(ii)Since fuzzy 𝑐-means clustering method is utilized for fuzzification in the proposed approach, the number of intervals, the interval length, and the degrees of memberships is not arbitrarily determined in the fuzzification phase, so the proposed approach is a systematic method.(iii)The methods in the literature that exploit fuzzy logic group relationship tables disregard the fuzzy set theory since only some elements of fuzzy sets are taken into account in the fuzzification phase. Thus, information loss occurs and it is expected to decrease the explanation power of the model. However, in the proposed method, particle swarm optimization method is used to calculate membership values in the fuzzy relationship matrix, so information loss and decrease in explanatory power of the model are prevented.(iv)In the proposed method, instead of using centroid method, a method which considers all membership degrees is employed. Therefore, information loss that occurs in the defuzzification stage is avoided.(v)Since the proposed approach provides the advantages given previously, it is expected that the suggested method has high forecasting accuracy.

The algorithm of the proposed method is given later step by step.

Algorithm 7. The proposed method’s algorithm.

Step 1. Fuzzify time series by using the fuzzy 𝑐-means clustering, and define fuzzy sets for absolute values of errors.
Let 𝑐 be the number of fuzzy set, such that 2≤𝑐≤𝑛. The fuzzy 𝑐-means clustering algorithm in which the number of fuzzy sets is 𝑐 is applied to the time series which consists of crisp values. After this application, the center of each fuzzy set is determined. Then, the degrees for each observation, which denote a degree of belonging to a fuzzy set for that observation, are calculated with respect to the obtained values of center of fuzzy sets. Finally, ordered fuzzy sets, 𝐿_𝑟(𝑟=1,2,…,𝑐), are obtained according to the ascending ordered centers, which are denoted by 𝑣_𝑟(𝑟=1,2,…,𝑐).
Fuzzy sets for absolute values of errors are defined in accordance with predetermined length of interval. Subintervals that contain absolute values of errors are generated. The number of clusters is taken as 𝑐 which equals the number of cluster of time series. Thus, the number of subintervals is 𝑐. The upper bound of the last subinterval is always open. While the first subinterval has the minimum absolute error, the subsequent subintervals have comparatively greater absolute errors. Generally, subintervals can be defined by 𝑢₁=[],𝑢0,𝑙𝑎₂=[𝑙𝑎,2×𝑙𝑎],…,𝑢_𝑐=[(𝑐−1)×𝑙𝑎,∞], (15) where 𝑙𝑎 represents the length of interval. In accordance with these subintervals, fuzzy sets for absolute errors are defined as follows: 𝐴₁=1𝑢₁+0.5𝑢₂,𝐴₂=0.5𝑢₁+1𝑢₂+0.5𝑢₃,𝐴_𝑐−1=0.5𝑢_c−2+1𝑢_c−1+0.5𝑢_c,𝐴_𝑐=0.5𝑢_c−1+1𝑢_c. (16)

Step 2. Determine the parameters of particle swarm optimization.
Some parameters of particle swarm optimization are possible intervals for social coefficient (𝑐_1𝑖,𝑐_1𝑓), cognitive coefficient (𝑐_2𝑖,𝑐_2𝑓), and inertia weight (𝑤₁,𝑤₂). The other parameters are swarm size (𝑝𝑛) and maximum iteration number (max𝑡). Evaluation function is RMSE criterion which is computing as in what follows: RMSE=∑^𝑛_𝑡=1𝑥_𝑡−𝑥_𝑡²𝑛, (17) where 𝑥_𝑡 is crisp time series, 𝑥_𝑡 is defuzzified forecasts, and 𝑛 is the number of forecasts.

Step 3. Generate a random initial population.
Elements of each particle are elements of 𝑅₁ and 𝑅₂ fuzzy relation matrices given in Definition 6. Since both 𝑅₁ and 𝑅₂ matrices consist of 𝑐² elements (𝑅_𝑘=[𝑟_𝑖𝑗𝑘],𝑖=1,2,…,𝑐,𝑗=1,2,…,𝑐,𝑘=1,2), the number of elements in a particle equals to 2𝑐². For the number of fuzzy sets 𝑐, representation of a particle is illustrated in Figure 1.
In this step, 𝑟_𝑖𝑗(𝑖,𝑗=1,2,…,𝑐) values of all particles are randomly generated from the interval [0,1] to generate a random initial population.

Figure 1 

The structure of a particle.

Step 4. Calculate the evaluation function values of all particles in the current swarm.
The method for calculating RMSE for any particle is given in Steps  4.1,  4.2, and 4.3.

Step  4.1. 𝑅₁ and 𝑅₂ fuzzy relation matrices are designed from positions of particles. Each 𝑐 position sets are the rows of 𝑅₁ and 𝑅₂ matrices.

Step  4.2. Fuzzy forecasts are calculated by using (14). After 𝑅₁ and 𝑅₂ relation matrices are estimated, fuzzy forecasts can be obtained by using fuzzy errors and fuzzy observations. When the first forecast is calculated, it is assumed that the error equals to zero since this assumption is also valid in conventional time series analysis.

For instance, some observations of a time series 𝑋_𝑡, corresponding membership values of these observations, and the center values (𝑣_𝑖,𝑖=1,2,3,4,5) for each fuzzy set are given in Table 1.

Table 1 

An example time series 𝑋_𝑡.

Let the relations matrices be as follows: 𝑅₁=10.5000001000000.500000010000,𝑅₂=.10.50000.501000.3000.0500000010000 (18)

Thus, fuzzy forecast for 1972 can be calculated as follows: 𝐹(1972)=𝐹(1971)∘𝑅₁∩𝜀(1971)∘𝑅₂==∘∩∘.0.70.50.2000.70.20.10010.5000001000000.50000001000010.500010.50000.501000.3000.0500000010000 (19)

After this procedure, fuzzy forecasts are defuzzified by using the rules given in what follows.(i)If the memberships of fuzzy forecast have only one maximum, then select the center of this cluster as the defuzzified forecasted value. (ii)If memberships of fuzzy forecast have two or more consecutive maximums, then select arithmetic mean of centers of corresponding clusters as the defuzzified forecasted value. (iii)Otherwise, standardize the fuzzy output, and use the center of the fuzzy sets as the forecasted value.

For example, defuzzified forecast for 1972 can be calculated as follows. Since the memberships of fuzzy forecast have only one maximum, defuzzified forecast value is taken as the center value (𝑣₁=13000) of the first fuzzy set which has the maximum membership value. This can be expressed by 𝑋₁₉₇₂=Def𝐹(1972)=13000. (20)

Then, absolute error value for 1972 is calculated given as follows: 𝐸₁₉₇₂=𝑋₁₉₇₂−𝑋₁₉₇₂=|13055−13000|=55. (21)

Absolute error value is fuzzified by using corresponding fuzzy set which has the maximum membership value for the interval which includes this absolute error value. For instance, for the time series given in Table 1, when the interval length of absolute error is taken as 50, 5 subintervals are as follows: 𝑢₁=[0,50],𝑢₂=[50,100],𝑢₂=[],𝑢100,150₂=[150,200],𝑢₅=[200,∞). (22)

Since the absolute error value 55 is in the second interval 𝑢₂, defuzzified absolute error is 𝐴₂. Thus, for 1972, the absolute error is obtained given as follows: 𝜀(1972)=[0.510.500]. (23) Then, the forecast for 1973 can be calculated as follows: 𝐹(1973)=𝐹(1972)∘𝑅₁∩𝜀(1972)∘𝑅₂. (24) In the same way, other forecasts for other years can be obtained.

Step  4.3. RMSE is computed using the formula given in (17).

Step 5. Update cognitive coefficient 𝑐₁, social coefficient 𝑐₂, and inertia parameter 𝑤 at each iteration by using formulas in (10).

Step 6. New positions of the particles are calculated by using the formulas given in (11) and (12).

Step 7. Repeat Steps 4 to 7 until maximum iteration bound (max𝑡) is reached.

5. Application of the Proposed Method

Firstly, the proposed method was applied to Alabama University enrollment data (1971–1992) that is well-known data in the literature. The enrollment data is presented in Table 2.

Table 2 

The enrollment data.

The algorithm of the proposed method was coded in MATLAB 7.9 version. In the application of the proposed method, we used seven fuzzy sets (𝑐=7) in Step 1 since the same number has been employed in other studies available in the literature (Cheng et al., 2008). The centers of clusters and the membership values which are obtained from FCM algorithm in Step 1 are given in Tables 3 and 4, respectively. In the fuzzification process of the absolute error, if the interval length is too small, aggregation will arise in the last fuzzy set. If it is too big, aggregation will occur in the first fuzzy set. Therefore, the length of interval for absolute error is picked as 300. Then, 7 fuzzy sets for absolute errors can be defined as follows: 𝑢₁=[0,300),𝑢₂=[300,600),𝑢₃=[),𝑢600,900₄=[900,1200),𝑢₅=[),𝑢1200,1500₆=[1500,1800),𝑢₇=[1800,∞). (25)

Table 3 

The centers of clusters obtained from FCM algorithm in Step 1.

Table 4 

The memberships of the observations obtained from FCM algorithm in Step 1.

In the literature, there are no general rules to determine parameter values of particle swarm optimization. Parameter values for this method have been generally specified intuitively due to the data in most of the applications. Therefore, the parameters of particle swarm optimization used in this study were intuitively determined like in other studies available in the literature (Ma et al., 2006). The parameters of the modified particle swarm optimization are determined according to trial and error method. The different values of the parameters were employed. The best parameters of the modified particle swarm optimization are determined as follows: (𝑐_1𝑖,𝑐_1𝑓)=(2,3),(𝑐_2𝑖,𝑐_2𝑓)=(2,3),(𝑤₁,𝑤₁)=(1,2),𝑝𝑛=30, and max𝑡=200. The optimal 𝑅₁ and 𝑅₂ matrices obtained from our methods are given later. The fuzzy and defuzzified forecasts of the proposed method are given in Table 5: 𝑅₁=𝑅1.00000.00001.00000.00001.00000.00000.45400.00000.75250.00000.04580.00000.00000.00000.33741.00001.00001.00001.00000.70001.00001.00000.00001.00001.00000.82830.00000.00000.42710.00001.00000.36270.00000.00001.00000.00001.00001.00001.00000.14431.00000.45790.00000.00001.00000.00001.00000.40020.0000₂,=1.00001.00000.14750.12730.58241.00001.00001.00000.23591.00000.64830.41690.26191.00000.70840.35261.00001.00000.00001.00001.00000.51221.00000.72120.00001.00001.00000.86760.00001.00001.00000.00000.11440.00000.00000.00000.00000.21080.00000.90190.00001.00000.31980.00000.14460.00000.00000.00000.6737. (26)(27)

Table 5 

The forecasts produced by the proposed methods.

The mean square error (MSE) value of the proposed method is given in Table 6. For comparison, the forecasting results of some other well-known fuzzy time series forecasting approaches are also presented in Table 6. When Table 6 was examined, it is seen that the proposed method gives the most accurate forecasts in terms of MSE criterion.

Table 6 

The obtained results for the enrollment data.

Secondly, to examine the forecasting performance of the proposed approach on test set, it was also applied to four different data sets. The first data set is index 100 in stocks and bonds exchange market of İstanbul (ISBEMI) whose daily observations are between May 20, 2008 and September 29, 2008. This data is called ISBEMI Set 1. The second data is index 100 in stocks and bonds exchange market of İstanbul (ISBEMI) whose daily observations are between October 3, 2008 and December 31, 2008. The second data is called ISBEMI Set 2. The third one is ISBEMI Set 3, and its daily observations are from October 01, 2009 to December 31 2009. ISBEMI Set 1, ISBEMI Set 2, and ISBEMI Set 3 have 95, 59, and 63 observations, respectively. For comparison, these time series were also forecasted using other well-known fuzzy time series method proposed by Song and Chissom [2], Chen [4], Chen [6], Huarng [5], Huarng and Yu [11], and Aladag et al. [29]. In the implementation, the last 15 of ESBEMI Set 1, 7 observations of ISBEMI Set 2, and 7 observations and 15 observations of ISBEMI Set 3 were used for test sets.

It is a well-known fact that there are no general rules to determine the parameters such as the length of interval or the number of fuzzy sets for some methods used in the implementation. Thus, related parameters were determined using trial and error method like in other studies available in the literature. To find the best length of interval for the methods introduced in [4], Chen [6], and Aladag et al. [29], the values between 200 and 2500 were examined with increment 100. The interval lengths in the methods given in Huarng [5] and Huarng and Yu [11] were determined when the methods were progressing because of the nature of these methods. The number of fuzzy sets was experienced between 5 and 25 for the proposed approach and Song and Chissom [2]. When the methods proposed by Chen [5] and Aladag et al. [29] were employed, models order 2, 3, 4, and 5 were examined since time series have daily observations. Finally, the number of neurons in hidden layer was changed from 1 to 10 when the approach proposed by Aladag et al. [29] was applied. After practicing, the forecasts obtained from the case where we have got the best result for the test data and the error criteria related to those forecasts are presented in Tables 7 and 8 for ISBEMI Sets 1 and 2, respectively.

Table 7 

The obtained results for ISBEMI Set 1.

Table 8 

The obtained results for ISBEMI Set 2.

For Table 7, the cases in which the superlative results were obtained are(i)the number of the fuzzy sets was 12 for the Song-Chissom method [2];(ii)the interval length was 1200 for the Chen method [4];(iii)the interval length was 800 for the Huarng distribution-based method [5];(iv)the interval length was 200 for the Huarng average-based method [5];(v)the ratio sample percentile was 0.5 for the Huarng and Yu ratio-based method [11];(vi)the number of the fuzzy sets was five for the Cheng et al. method [19];(vii)the number of the fuzzy sets was 11 and the number of the neurons in the hidden layer was five for Yolcu et al. [32];(viii)the number of the fuzzy sets was 5 in the proposed method.

For Table 8, the cases in which the superlative results were obtained are(i)when the number of fuzzy sets is 12, for the method of Song and Chissom [2];(ii)when the length of interval is 1900 for the method of [4];(iii)when the length of interval is 2200 and the model order is 2 for the method of Chen [6];(iv)when the length of interval is 800 for the method of Huarng distribution based [5];(v)when the length of interval is 200 for the method of Huarng average based [5];(vi)when the ratio sample percentile is 0.5 for the method of Huarng and Yu ratio based [11];(vii)when the length of interval is 1000, the model order is 2, and number of neurons in hidden layer is 7 for the method proposed by Aladag et al. [29];(viii)when the number of fuzzy sets is 5 in the proposed method.

The best cases for the results given in Table 9 resulted when(i)Song-Chissom method [2] was applied with nine fuzzy sets;(ii)interval length was 1300 for the Chen method [4];(iii)interval length was 800 for the Huarng distribution-based method [5];(iv)interval length was 200 for the Huarng average-based method [5];(v)ratio sample percentile was 0.50 for the Huarng and Yu ratio-based method [11];(vi)the number of the fuzzy sets was 13 and the number of the neurons in the hidden layer was seven for Yolcu et al. [32];(vii)when the number of fuzzy sets is 5 in the proposed method.

Table 9 

The obtained results (last 7 observations are test set) for ISBEMI Set 3.

The best cases for the results given in Table 10 resulted when(i)Song-Chissom [2] method was applied with nine fuzzy sets;(ii)interval length was 1500 for the Chen method [4];(iii)interval length was 800 for the Huarng distribution-based method [5];(iv)interval length was 200 for the Huarng average-based method [5];(v)ratio sample percentile was 0.50 for the Huarng and Yu ratio-based method [11];(vi)the number of the fuzzy sets was seven and the number of the neurons in the hidden layer was three for Yolcu et al. [32];(vii)the number of the fuzzy sets was five in the proposed method.

Table 10 

The obtained results (last 15 observations are test set) for ISBEMI Set 3.

When Tables 6–10 are examined, it is obvious that the most accurate forecasts are obtained for all data sets when the proposed forecasting approach is used.

6. Conclusion and Discussion

In this study, an ARMA type fuzzy time series forecasting model is firstly introduced in the literature. In addition, to analyze this model, a novel approach based on fuzzy 𝑐-means and particle swarm optimization methods is suggested in this study. To show the forecasting performance of the proposed forecasting approach, the enrollments of Alabama University are applied to a well-known data and obtained forecasting results are compared to those produced by other fuzzy time series forecasting methods available in the literature. In addition, to evaluate the forecasting performance of the proposed forecasting approach on test set, two different time series are forecasted by utilizing the proposed method and other approaches available in the literature. As a result of the comparison, it is observed that the proposed approach gives the most accurate forecasts.

In the literature, the proposed method is the first method based on a fuzzy time series forecasting model which contains both AR and MA terms. In the linear stochastic ARMA models, the AR(∞) model is equal to MA(1) model. But this result is only valid for linear ARMA models. There is no such kind of result for fuzzy time series model, artificial neural network model, and nonlinear stochastic models. Because of this, adding MA terms to fuzzy time series forecasting model is really important. Also, in the proposed method, particle swarm optimization method is used to calculate membership values in the fuzzy relationship matrix, so information loss and decrease in explanatory power of the model are prevented. Besides, since fuzzy 𝑐-means clustering method is employed in the fuzzification process, the number of intervals, the interval length, and the degrees of memberships are not arbitrarily determined, so the proposed approach works in a systematic way. In addition, information loss which occurs in the defuzzification stage is avoided by using a method which takes all membership degrees into account instead of using centroid method. Therefore, the proposed method has high forecasting accuracy.

Moreover, long-range dependence was determined for all time series by using 𝑅/𝑆 tests, and results are given in the appendix. All of the time series are short-range dependent. In the future studies, the performance of the proposed method will be able to research for long-range time series. But in the proposed fuzzy time series method, any linear model was not used. Fuzzy time series methods can be used to forecast nonlinear time series. In the fuzzy time series methods, the relations are employed instead of functions.

The proposed forecasting approach is suggested to forecast a first-order model including both AR and MA terms. In future studies, high-order models can be defined. Then, it is possible to extend the proposed approach to forecast high order models.

Appendix

𝑅/𝑆 test was applied to five time series to determine long-range dependence. The test results are given in Table 11. As a result of hypotheses tests, all data sets are short-range dependent. The tests were applied in S-PLUS package program.

Table 11 

Results of 𝑅/𝑆 test.

Acknowledgments

The authors would like to thank the reviewers for their helpful comments and opinions. This work was supported by “The Scientific and Technological Research Council of Turkey (TUBITAK),” Turkey, under Project no. 210T150.