Fuzzy Time Series: An Application to Tourism Demand Forecasting

: Problem statement: Forecasting is very important in many types of organizations since predictions of future events must be incorporated into the decision-making process. In the case of tourism demand, better forecast would help directors and investors make operational, tactical and strategic decisions. Besides that, government bodies need accurate tourism demand forecasts to plan required tourism infrastructures, such as accommodation site planning and transportation development, among other needs. There are many types of forecasting methods. Generally, time series forecasting can be divided into classical method and modern methods. Recent studies show that the newer and more advanced forecasting techniques tend to result in improved forecast accuracy, but no clear evidence shows that any one model can consistently outperform other models in the forecasting competition. Approach: In this study, the performance of forecasting between classical methods (Box-Jenkins methods Seasonal Auto-Regressive Integrated Moving Average (SARIMA), Holt Winters and time series regression) and modern methods (fuzzy time series) has been compared by using data of tourist arrivals to Bali and Soekarno-Hatta gate in Indonesia as case study. Results: The empirical results show that modern methods give more accurate forecasts compare to classical methods. Chen’s fuzzy time series method outperforms all the classical methods and others more advance fuzzy time series methods. We also found that the performance of fuzzy time series methods can be improve by using transformed data. Conclusion: It is found that the best method to forecast the tourist arrivals to Bali and Soekarno-Hatta was to be the FTS i.e., method after using data transformation. Although this method known to be the simplest or conventional methods of FTS, yet this result should not be odd since several previous studies also have shown that simple method could outperform more advance or complicated methods.


INTRODUCTION
Forecasting is very important in many types of organizations since predictions of future events must be incorporated into the decision-making process. In the case of tourism demand, better forecast would help directors and investors make operational, tactical and strategic decisions. Besides that, government bodies need accurate tourism demand forecasts to plan required tourism infrastructures, such as accommodation site planning and transportation development, among other needs. There are many types of forecasting methods. Generally, in time series we can divide forecasting method into classical or traditional method and modern methods. Although recent studies show that the newer and more advanced forecasting techniques tend to result in improved forecast accuracy under certain circumstances, no clear-cut evidence shows that any one model can consistently outperform other models in the forecasting competition (Song, 2008).
The time series forecasting methods have found applications in very wide areas including finance and business, computer science, all branches of engineering, medicine, physics, chemistry and many interdisciplinary fields. Conventionally, researchers have employed traditional methods of time series analysis, modeling and forecasting. Some of mainly been used that will discuss in this study are Box-Jenkins methods Seasonal Auto-Regressive Integrated Moving Average (SARIMA), Holt Winters and time series regression. The conventional time series modeling methods have served the scientific community for a long time. However, they provide only reasonable accuracy and suffer from the assumptions of stationarity and linearity. Due to these constrains, comes the idea of alternative solution that is fuzzy time series. In this study, the performance of forecasting between classical (Box-Jenkins methods Seasonal Auto-Regressive Integrated Moving Average (SARIMA), Holt Winters and time series regression) and fuzzy time series has been compared. Song and Chissom (1993a;1993b) first introduced the definitions of fuzzy time series and developed their model by using fuzzy relation equations and approximate reasoning. Since that, fuzzy time series has gains much attention from researchers in many fields and the methods have been developed rapidly. In forecasting time series, Chen (1996) proposed a first order fuzzy time series used simplified arithmetic operations and fuzzy logical relationship groups to forecast the enrollments of the University of Alabama. Then, Chen (2002) developed a high-order fuzzy time series model by extending Chen's first-order model (Chen, 1996). Lee and Chou (2004) improved the forecasting accuracy of Chen's model (Chen, 1996) by properly defining the number of linguistic variables. In order to overcome recurrence and weighting problems in fuzzy time series forecasting, Yu (2005) developed the weighted fuzzy time series models. Recently, in tourism demand forecasting Tsaur and Kuo (2011) developed adaptive fuzzy time series model to forecast Taiwan's tourism demand. Suhartono and Lee (2011) proposed a hybrid approach based on winter's model and weighted fuzzy time series in order to analyze trend and seasonal fuzzy time series for tourism data to Bali.
In this study we used data of tourist arrivals to Bali and Soekarno-Hatta gates in Indonesia as case-study. The data were taken from the Indonesia Central Bureau of Statistics. All the dataset contain monthly data from January 1989 to December 1997. We only consider the data until 1997 to anticipate extreme data (Bali bombing). For the estimation (in-sample) purpose, data are taken from January 1989 to December 1996. Meanwhile, data from January 1997 to December 1997 are considered for the testing or evaluation (out-sample) purpose.

MATERIALS AND METHODS
The SARIMA model: The Box-Jenkins approach to modelling autoregressive integrated moving average (ARIMA) processes involved an iterative three-stage process of model selection or identification, parameter estimation and model checking.
Since the tourist arrivals data that we used in were measured at regular calendar intervals within a year, it may exhibit periodic behaviour. Hence, the general Box-Jenkins model which allocates seasonality with P seasonal autoregressive terms, D seasonal differences and Q seasonal moving average terms (Box et al., (1994) is given as follows: The above equation also can be written as follows (Bowerman et al., 2005): The purpose of dummy variable is to ensure that an appropriate seasonal parameter is included in the regression model in each time period. This dummy variable model assumes that the seasonal component is unchanging from year to year.

Time series regression (multiplicative model):
Regression models involving trigonometric terms can be used to forecast time series exhibiting increasing seasonal variation (Bowerman and O'Connell, 1993): This model also has a linear trend assumption and it was altered to handle trends in our case study. Hence, the following multiplicative model was obtained: Holt Winters: Method proposed by Winters in 1960 is one of the most famous forecasting techniques for seasonal time series. It has two types of method, additive seasonality and multiplicative seasonality. In this study we use the multiplicative Holt Winters method which is more often applied in software package. The magnitude of the seasonal variation increase when the mean level of the time series increase, hence the seasonality is multiplicative. The three smoothing equations of the multiplicative Holt Winters method given by the following equations: Where: L t = Level at time t T t = Trend at time t S t = Seasonal component at time t y t = Fitted value or one-period-ahead forecast at time t y = Weight for the trend a = Weight for the level y t = Data value at time t δ = Weight for the seasonal component P = Seasonal period The parameters α, β, γ should lie in the interval (0, 1). We use the default weight in Minitab, which the weight for this level, trend and seasonal component are 0.2. Song and Chissom (1993b) first introduced the definitions of fuzzy time series and developed their model by using fuzzy relation equations and approximate reasoning. General definitions of fuzzy time series are given as follows.

Fuzzy time series:
Let U be the universe of discourse, where U = {u 1 , u 2 ,…u b }. A fuzzy set A j of U is defined as A j = f Aj (u 1 )/u 1 +f Aj (u 2 )/u 2 +…+ f Aj (u b )/u b , where F Aj is the membership function of the fuzzy set A j ; f Aj :U→[0, 1] u a is a generic element of fuzzy set A j ; f Aj (u a ) is the degree of belongingness of u a to A j ; f Aj (u a ) ∈[0, 1] and 1≤a≤b.
Definition 1: Fuzzy time series. Let y (t) (t = …, 0, 1, 2…) a subset of real numbers R, be the universe of discourse by which fuzzy sets f i (t) are defined. If F(t) is a collection of f 1 (t), f 2 (t),… then F(t) is called a fuzzy time series defined on Y(t).

Definition 2:
If there exists a fuzzy relationship R (t-1, t), such that F (t) = F (t-1) R (t-1, t), where ° is an arithmetic operator, then F (t) is said to be caused by F (t-1). The relationship between F (t) and F (t-1) can be denoted by: F (t-1) →F (t).
Definition 4: Suppose F (t-1) = A j and F (t) = A j , a fuzzy logical relationship can be defined as: where, A i and A j are called the Left-Hand Side (LHS) and Right-Hand Side (RHS) of the fuzzy logical relationship, respectively. The definition of the first order seasonal fuzzy time series model for forecasting proposed by Song and Chissom (1993b) is given as follows.
Definition 5: Let F(t) be a fuzzy time series. Assume there exists seasonality in {F(t)}, first order seasonal fuzzy time series forecasting model: where, m denotes the period.
The high order fuzzy time series model proposed by Chen (2002) is given as follows: Definition 6: Let F(t) be a fuzzy time series. If F(t) is caused by F (t-1), F (t-2,… and F (t-n) then this fuzzy logical relationship is represented by: And it is called the n th order fuzzy time series forecasting model.
Initially, the repeated FLRs were simply ignored when fuzzy relationships were established. In many previous studies, each FLR was treated as if it was of equal importance, which may not have reflected the real world situation (Chen, 1996;Huarng, 2001;Song and Chissom, 1993a;1993b;Yu, 2005). In this scenario, the occurrences of the same FLRs are regarded as if there were only one occurrence. In other words, the recent identical FLRs are simply ignored. To explain this, suppose there are FLRs in chronological order that have the same LHS.
A 1 , as follows Eq. 1: Following Chen (1996), these FLRs in Eq. 1 are used to establish an FLRG as: The ignoring of recurrence, however, is questionable. Yu (2005) argued that the occurrence of a particular FLR represents the number of its appearances in the past. For instance, in Eq. 1, A 1 →A 1 appears three times, both A 1 →A 2 and A 1 →A 2 only once. The recurrence can be used to indicate how the FLR may appear in the future.
Later, Yu (2005) proposed the chronological weights to deal with recurrent fuzzy relationships and their importance. To illustrate it, suppose there are FLRs in chronological order as in Eq. 1 and then the weights are as follows: As a result, the most recent FLR (t = 5) is assigned the highest weight of 5, which means that the probability of its appearance in the near future is higher than in the case of the others. On the other hand, the most aged FLR (t = 1) is assigned the lowest weight of 1, which means that the probability of its appearance in the near future is lower than in the case of the others.
Recently, Cheng et al. (2008) proposed the weights focused on the probability of its appearance and their importance of chronological FLR for the same recent identical FLRs. To explain it, suppose there are FLRs in chronological order as in Eq. 1 and then the weights are as follows: In this study, we applied the FTS according to the following procedures: Step 1: The data are not stationary, hence data preprocessing has been carried out by taking transformation according to Koehler Eq. 2: where, λ is the coefficient from Box-Cox Transformation.
Step 2: Fuzzy relationship was determined according to SARIMA model for the data set. This procedure also has been done by Faraway and Chatfield (1998) in order to select neural network input variable. Due to limitation of considered methods, we have to ignore the lag from MA and SMA. For instance SARIMA model for Bali is (0, 1, 1) (0, 1, 1) 12 (Hence the fuzzy relationship can be denoted by: Step 3: In order to select input and fuzzy time series order, firstly we try to input all the three input from fuzzy relationship that obtained in step 2.
To experiment the selection of input, we try all the possible combination of two input from the three inputs and single input as well. So all possible input are Lag 1, 12, 13; Lag 1, 12; Lag 1, 13; Lag 12, 13; Lag 1; Lag 12; and Lag 13.
Step 4: The optimum length of intervals was calculated following average-based length by Huarng (2001).
Step 6: Forecast data were transformed back to original scale data and the forecast accuracy were calculated.  For in sample, n is the number of observations (degree of freedom in case of SARIMA). Meanwhile, for out sample, n is the number of forecast data which is 12 in this study.

SARIMA:
In this study we used MINITAB version 14 to analyze SARIMA model. In model identification stage, firstly we used time series plot to see briefly whether the data have seasonal and trend patterns. Figure 1 and 2 show clearly that both data sets have trend patterns; hence the assumption of stationary condition in mean is not satisfied. Yet, we validate this assumption by using Autocorrelation Function (ACF) and Partial Autocorrelation (PACF) plot and from the results ACF plot for both data only dies down until first differencing both in non-seasonal and seasonal (S = 12). From box-cox plot, we found out that both data are not stationary in variance. Data transformation was made according to Johnson and Wichern (2007). After the data preparation, we plot again ACF and PACF in order to identify the model and took a few possible models. The best model was chosen among these competitive models base on the smallest RMSE value. The model parameters were estimated by using least squares estimation. Models with insignificant parameters (exclude constant) were eliminated. Remaining models then proceed to the diagnostic checking step to see whether the models adequate by using Ljung-Box statistics. According to RMSE value, the best SARIMA models for Bali is SARIMA (0,1,1) (0,1,1,) 12 , meanwhile for Soekarno-Hatta is SARIMA (1,1,0) (1,1,1,) 12 .Therefore, the models to forecast tourist arrivals to Bali and Soekarno-Hatta after taking the parameter estimation from MINITAB output are following the below equations respectively: Fuzzy time series: All the three methods of FTS were implemented by using MatlabR2008a. The selection of the best FTS input is according RMSE. We choose only lag 13 for Bali and lag 2, 3, 13, 14, 15 and 25 for Soekarno-Hatta. In case of Chen's method for Bali, since only one lag was chosen, hence it follows Chen's first order method Chen (1996). The number of intervals for Bali and Soekarno-Hatta are 26 and 30 respectively (calculated following average-based length procedure). Figure 1 and 2 show time series plot for tourist arrivals to Bali and Soekarno-Hatta respectively before data transformation. Meanwhile, Figure 3 and 4 show after data transformation following Eq. 2 for tourist arrivals to Bali and Soekarno-Hatta respectively. The Comparison of forecast performance between three F.T.S methods before and after data transformation is shown in Table 1. Before the data transformation, Cheng's methods outperformed all other methods, but after transformation Chen's methods is the best according to MAPE, MAD and RMSE. There is exception for Bali before data transformation, where Chen's method is the best according to RMSE. Comparison of forecasting performance: The forecasting of tourist arrivals to Bali and Soekarno-Hatta gate in testing period is done using four different approaches (total seven methods). From the results in Table 2, it is clearly shows that FTS Chen's and FTS Yu's method outperform all the classical methods and FTS Chen's method gives the most accurate forecast according to MAPE, MAD and RMSE. It seem that all the three accuracy measurement give inconsistent ranking, except for the first ranking which Chen's methods has been chosen by all of the measurement.

DISCUSSION
Compare with Fig. 1 and 2, it seem that the transformation produced series with more constant variations. The changing rank before and after the data transformation might due to different assign weight for FLR by both methods, since Chen's method gives the same weight and Cheng's method gives different weight according to the number of times each FLR occur.    Despite the changing rank, when we compare the magnitude value of error, the data transformation still improved the forecast accuracy for the best method of both data sets. For tourist arrivals to Bali, Holt Winter's has lowest value of RMSE after Chen's method. Holt Winter's method which is exponential smoothing model is basically weight average of previous observations. Hence it can predict well for series that repeat their pattern and scale like Bali

CONCLUSION
It is found that the best method to forecast the tourist arrivals to Bali and Soekarno-Hatta was to be the FTS i.e., Chen (1996;2002) method after using data transformation. Although this method known to be the simplest or conventional methods of FTS, yet this result should not be odd since several previous studies also have shown that simple method could outperform more advance or complicated methods (Makridakis et al. (1982;; Spyros Makridakis and Hibon (2001). According to Spyros Makridakis and Hibon (2001) many simple methods, such as a random-walk model, for example, offer adaptability to structural change, in that the model immediately adapts to the latest level of the series.  (1997) Generally, we can conclude that FTS (especially Chen's method) is good to predict fluctuating series such as tourist arrivals. However, future research could be done by using FTS that can consider Moving Average (MA) terms in input lag in order to improve the forecast accuracy. More advance accuracy measurement such as statistical test also can be applied to evaluate the forecast accuracy among competition models. Besides, turning points and directional change errors also will be useful since they can give better information on tourism growth cycles instead of just forecast error magnitude. When refer to Fig. 5 and 6, all the consider methods perform better for Bali compare to Soekarno-Hatta. This is due to the different pattern of both data. Bali continue the pattern of training data in testing data, in contrast this is not happen to Soekarno-Hatta which dropped dramatically in May 1997 and then increase sharply until July 1997. This shows that tourist arrivals to Bali has less affected by economic crises in Asia compare to Soekarno-Hatta.