A new explainable robust high-order intuitionistic fuzzy time-series method

Fuzzy time series, based on type-1 fuzzy sets, continue to have a wide range of use in the literature. These methods use only membership values to determine the fuzzy relations. However, intuitionistic fuzzy time series models use both membership values and non-membership values. So it can be considered that the use of intuitionistic fuzzy time forecasting models will be able to increase the forecasting performance because the intuitionistic fuzzy sets have more information than fuzzy sets. Therefore, intuitionistic fuzzy time series models have started to employ for forecasting the real-life series in the fuzzy time series literature. A novel, explainable, robust high-order intuitionistic fuzzy time series forecasting method is proposed based on a newly defined model. In the proposed method, the intuitionistic fuzzy c-means algorithm is used for the fuzzification of observations, and a robust regression method employed for determining fuzzy relations. With the use of robust regression in determining the fuzzy relationships, all inputs of the proposed method can be explainable and they can be tested and commented on statistically. Applications of this study are made by using energy data of Primary Energy Consumption between the years 1965 and 2016 for 23 countries in the region of Europe-Eurasia. The forecasting performance of the proposed method is compared with the performance of some selected benchmarks, and the obtained results are discussed.


Introduction
Forecasting techniques continue to be largely used in the time series analysis at present. One of the most popular forecasting techniques is the autoregressive integrated moving average (ARIMA) (Kadılar 2005). The approach of Box-Jenkins (1976) is also frequently used as a basic approach in time series literature to obtain the best ARIMA model. In addition, exponential smoothing methods have also been used in the literature for time series forecasting. Indeed, exponential smoothing models are special cases of ARIMA models. In classical time series, exponential smoothing methods are generally described under three headings as the simple exponential smoothing, Holt's exponential smoothing and Winter's exponential smoothing (Kadılar 2005). However, there are also new exponential smoothing methods proposed in the literature. One of these methods is ATA exponential smoothing method proposed by Yapar (2016). ATA exponential smoothing method (Yapar 2016;Yapar et al. 2017a, b) is an updated and effective time series analysis method proposed as an alternative to exponential smoothing methods in the literature.
Conventional forecasting techniques need assumptions like normality, stability and linearity because these methods are based on probability theory. However, many time series in real life do not satisfy any assumptions, and so these time series need to be solved via using some nonprobabilistic methods in the literature. In the literature, non-probabilistic time series forecasting methods are Communicated by Mu-Yen Chen. & Erol Egrioglu erol.egrioglu@giresun.edu.tr; e.egrioglu@lancaster.ac.uk generally separated into three parts as fuzzy time series, artificial neural networks and other computational methods. Additionally, these non-probabilistic methods could be used together. For instance, the analysis of a fuzzy time series method includes an artificial neural network or a computational method. In this regard, the content of this study is generally about a new non-probabilistic method. The first fuzzy time series approach was proposed by Song and Chissom (1993a). In Song and Chissom (1993a), the analysis of a fuzzy time series consists generally of three phases. These phases are fuzzification, the determination of fuzzy relations and defuzzification, respectively. In the fuzzification phase, dividing the universe of discourse proposed by Song and Chissom (1993a, b), was used largely in literature. Additionally, there are some studies about being selected for the optimal interval length in the literature because determining the interval lengths affects the forecasting performance significantly in the fuzzification phase. In the literature, the first study of these studies was proposed by Huarng (2001) to determine the optimal interval length. Huarng (2001) proposed two different approaches based on distribution and average. Huarng and Yu (2006) also proposed a ratio-based method to determine the optimal interval length. The other studies to determine the optimal intervals were given by Yolcu et al. (2009), Egrioglu et al. (2010) and Egrioglu et al. (2011). In fuzzy time series literature, the centralization technique is mainly used in the defuzzification phase that is the final phase to solve the fuzzy time series. The determining of fuzzy relations is very important in the fuzzy time series methods. Therefore, this stage is effective to improve the forecasting performance of any fuzzy time series method. Song and Chissom (1993b) used the complex matrix operations to determine fuzzy relations as being based on the fuzzy set theory of Zadeh (1965). Then, Chen (1996) proposed a new method based on fuzzy relation tables instead of the complex matrix operations of Song and Chissom (1993b). Thereby, many fuzzy time series methods based on fuzzy relation tables have been proposed so far. In fuzzy time series literature, some methods that use fuzzy relation tables to determine fuzzy relations were given by Chen (2002), , Chen and Tanuwijaya (2011), Kocak (2013), Cheng et al. (2016) and Kocak (2017). Artificial neural networks (ANNs) are intensely used in fuzzy time series. Aladag et al. (2009) used the ANNs to determine fuzzy relations for the high-order fuzzy time series model in the fuzzy time series literature. Egrioglu et al. (2009a) also used the ANNs for the high-order multivariate time series method. Yu and Huarng (2010) proposed a new method based on the artificial neural network to define fuzzy relations. In fuzzy time series literature, some methods that use ANNs to determine fuzzy relations were given by Egrioglu et al. (2009b, Aladag (2013), Kocak (2015), Bas et al. (2015), Kocak (2013) and Bas et al. (2018). Also, many fuzzy time series methods have been proposed based on particle swarm optimization (PSO). In the fuzzy time series literature, Kuo et al. (2009) developed a new method named hybrid particle swarm optimization (MPSO) for solving the historical data of enrollments of the University of Alabama. Similarly, Hsu et al. (2010) used the modified PSO method and proposed a new method named modified turbulent particle swarm optimization (MTPSO) for forecasting temperature. Some methods that use PSO were also given by Kuo et al. (2010), Park et al. (2010, Huang et al. (2011), Aladag et al. (2012, Chen and Kao (2013), Singh and Borah (2014) and Cagcag Yolcu and Lam (2017). Besides, some studies use fuzzy c-means (FCM) to obtain membership values. Egrioglu et al. (2013) and Cagcag Yolcu and Alpaslan (2018) proposed different models based on FCM and PSO. There are also non-probabilistic and non-parametric time series forecasting methods based on the methods of machine learning, fuzzy inference system, artificial neural network, artificial intelligence, besides fuzzy time series methods. Some of these studies were given by Yolcu et al. (2013), Aladag et al. (2014), Zhou et al. (2016), Egrioglu et al. (2017) and Tak (2021). Yolcu et al. (2013) proposed a new artificial neural network based on PSO for solving linear and nonlinear time series. Aladag et al. (2014) proposed a new artificial neural network time series forecasting method based on the median neuron model. Zhou et al. (2016) predicted some financial time series by using a dendritic neuron model. Tak (2018) proposed a meta fuzzy functions approach for time series forecasting. Tak et al. (2018) proposed a recurrent type-1 fuzzy functions approach for time series forecasting. Egrioglu et al. (2019a) proposed a median-Pi ANN model for time series analysis. Tak (2020a) proposed a type-1 possibilistic fuzzy forecasting functions method for time series forecasting. Fuzzy time series uses only membership values for solving fuzzy time series. However, intuitionistic fuzzy time series based on intuitionistic fuzzy sets use both membership and non-membership values. Therefore, more information is used for the forecasting process. The first article about intuitionistic fuzzy time series (IFTS) was proposed by Zheng et al. (2013a). For this reason, the subject of IFTS can be thought a new and an updated field for time series literature. Fundamental definitions of IFTS were initially given by Zheng et al. (2013a). Furthermore, Zheng et al. (2013b) used the intuitionistic fuzzy c-means algorithm proposed by Chaira (2011) to determine membership and non-membership values. Lei et al. (2016) proposed a high-order intuitionistic fuzzy time series model by using an unequal partition of the universe of discourse based on the fuzzy clustering algorithm. Kumar and Gangwar (2016) proposed a method to determine intuitionistic fuzzy relations. Fan et al. (2016) proposed a long-term intuitionistic fuzzy time series forecasting model. Egrioglu et al. (2019b) proposed an intuitionistic fuzzy time series forecasting method. Also, some other studies about intuitionistic fuzzy time series were proposed by Zheng et al. (2014), Wang et al. (2014), Joshi et al. (2016), Hu et al. (2017), Fan et al. (2017), Abhishekh et al. (2018 and Tak (2020b Let IF t be an intuitionistic fuzzy time series. A 1 ; A 2 ; . . .; A c are intuitionistic fuzzy sets on a universal set. l A j t ð Þ, m A j t ð Þ are membership and non-membership values of tth observation to jth intuitionistic fuzzy set. The high-order single-variable intuitionistic fuzzy time series model is given in Eq. (1).
In Eq. (1), G is a linear or nonlinear function, e t is an error term with zero mean, l A 1 t À 1 ð Þ; . . .; l A c t À p ð Þ and m A 2 t À 1 ð Þ; . . .; m A c t À p ð Þ are lagged membership values and non-membership values, respectively, obtained from X tÀ1 ; . . .; X tÀp . Equation (1) expresses that any time series are affected by lagged membership values and nonmembership values.
Definition 2 (New Definition: High-Order Single-Variable Intuitionistic Fuzzy Time Series Forecasting Model based on the principal components analysis and the robust regression analysis) Let IF t be an intuitionistic fuzzy time series. A 1 ; A 2 ; . . .; A c are intuitionistic fuzzy sets on a universal set. l A j t ð Þ, m A j t ð Þ are membership and non-membership values of tth observation to jth intuitionistic fuzzy set for j ¼ 1; 2; . . .; c. The high-order single-variable intuitionistic fuzzy time series model based on the principal components analysis is given in Eq. (2).
In Eq. (3), p is the number of lagged variables, l A 1 t À 1 ð Þ; . . .; l A c t À p ð Þ and m A 1 t À 1 ð Þ; . . .; m A c t À p ð Þ are lagged membership values and non-membership values obtained from X tÀ1 ; . . .; X tÀp , q is the number of the important principal component for q cp, f i is ith linear function obtained by using the correlation matrix of the principal component analysis,z i is the ith z-score calculated via the principal component analysis, b 0 ; b 1 ; b 2 . . .b q are regression coefficients obtained via the robust regression method and e t is the error term with zero means.
3 Intuitionistic fuzzy C-means algorithm Chaira (2011) proposed a clustering algorithm (IFCM) for intuitionistic fuzzy sets. The algorithms have been commonly used in intuitionistic fuzzy time-series studies (Zheng 2013a, b;Egrioglu et al. 2019b) in the literature. The algorithm of intuitionistic fuzzy c-means clustering was given in Algorithm 1.

Algorithm 1
The intuitionistic fuzzy c-means clustering algorithm.
In Eq. (7), x k k ¼ 1; 2; . . .; n ð Þare observations of the time series the and f is the fuzziness index.
In Eq. (10), u Ã ik and p ik are intuitionistic membership values and hesitation degrees obtained in Step 5.
Step 7. Check the stopping criteria given in Eq. (11). If the condition is satisfied, stop the algorithm, otherwise make U old ¼ U new and go to Step 3. e is a small positive number and k:k 2 is the L 2 norm.
In Eq. (11), e is a small positive number and k:k 2 is the L 2 norm.

The robust regression
The robust regression is commonly used instead of the least-square method in the literature because robust estimators are not affected by outliers too much. Robust regression techniques can be given as the L1 technique of Edgeworth (1887), M-Estimates of Mallows (1975), and S-Estimates of Rousseeuw (1984). Robust regression algorithms based on these fundamental robust techniques have been proposed in the literature. One of these algorithms is the iteratively reweighted least-squares technique (James et al. 1988). In this study, robust regression estimates were used iteratively reweighted the least-squares technique in defining the fuzzy relation phase. Accordingly, the algorithm of the iteratively reweighted least squares is given in Algorithm 2.
Step 1. Calculate initial regression coefficient estimates b b of b via the least-square method for linear regression model given in Eq. (12) and calculate standard deviation estimate b r of r by using Eq. (13).
In Eq. (12), n is the observation number, p is the independent variable number, Y is the dependent variable vector, and X is the independent variable matrix that values of its first column are 1, e is an error parameter vector, e i is an estimate of e i for ith observation obtained via the least and 0.6745 is the constant used to make estimates for normal distribution.
In Eq. (14), y i is ith observation value of the dependent variable and x 0 i is a row vector that is the transpose of the first column of X matrix.
In Eq. (15), b r is standard deviation estimates and 4.685 is the tuning constant for the bisquare distribution.
Step 6. Check the stopping criteria given in Eq. (19). If the condition is satisfied, stop the algorithm, otherwise In Eq. (19), e is a small positive number and : j j is the absolute value.

The proposed method
In this study, a new intuitionistic, high-order fuzzy time series forecasting method was proposed. The proposed method generally uses Eq. (3) given in Definition 1 and completely uses Eqs. (4)-(5) given in Definition 2. Thereby, both membership and non-membership are used in the intuitionistic fuzzy time series approach, unlike the fuzzy time series approach that uses only non-membership values. In the proposed method, the intuitionistic fuzzy c-means clustering method is used in the fuzzification phase. Furthermore, the principal component analysis is used to eliminate the multicollinearity problem and reduce dimension because there are many inputs before the phase of determining fuzzy relationships of the proposed method. The proposed method is based on the robust regression. For this reason, the iteratively reweighted least-squares technique is a robust regression method that is used to define functional relationships in the stage of determining the fuzzy relationships. Also, the proposed method does not need the defuzzification step because the outputs of the proposed method are real forecasts. The flowchart and algorithm of the proposed method is given in Fig. 1 and Algorithm 3, respectively.
Algorithm 3 Algorithm of the proposed method.
Step 1. Determine X t , X test t ; X train t , n; n test , n train ; c and m ð1; 2; 3; . . .Þ are data and parameters used in the application of the proposed method. Here, X t is the time series, X test t is the test set ofX t , X train t is the training test ofX t , n is the length of X t , n test (the last part of %5, %10 or %15 of n) is the length of X test t , n train is the length of X train t , c (2,3,…) is the number of fuzzy sets and m (1,2,3,…) is the order of the model.
Step 2. Calculate membership and non-membership for the training time series (X train t Þ by using the intuitionistic fuzzy c-means (IFCM) clustering method via Algorithm 1 given in Sect. 3.
Step 3. Apply lag process for U and V separately. And so, obtain input matrix ( X input ) that includes lagged variables of membership and non-membership under Eq. (5) given in Definition 2.
For example, X input matrix obtained by using U and V in Eqs. (20) and (21) is given in Eq. (24).
Step 4. Principal component analysis (PCA) is applied to X input matrix obtained in Step 3. And so, obtain the matrix of principal component coefficients (PCC), the matrix of principal component (PC), the matrix of important principal component (IPC) and the new input matrix ðP input Þ; respectively.
Firstly, elements of X input matrix are transformed to the standard normal matrix named Z input that means and standard deviations of each column of Z input are 0 and 1. Then, PCA is applied to Z input . And so, PCC and PC matrix are calculated. Also, important principal components are selected from principal components by taking into consideration the criterion that the explained variance ratio needs to be at least 0.85. In this way, IPC matrix is obtained. Lastly, the new input matrix named P input is obtained by putting a vector that consists of 1 value to the first column of IPC.
For example, PCC, PC, IPC and P input matrices obtained by applying PCA to X input in Eq. (25) are given below.
In Eq. (25), '' T '' indicates transpose operation. The explained variance ratio of two important principal components is %97.86. So, it can be said that the aims of dimension reduction and obtaining orthogonal inputs are accomplished via PCA.
Step 5. Define fuzzy relations by using robust Regression.
Robust regression coefficient estimates ( b b) is calculated via iteratively reweighted least-squares algorithm given in Algorithm 2 by taking X ¼ P input and Y ¼ X train For example, b b obtained via Algorithm 2 by taking X ¼ P input in (29) and Y ¼ X train t ¼ À; À; 11; 12; 16; 13 ½ is given below.
Step 6. Apply determining membership values for time series (X t Þ to obtain new memberships (U new Þ and new non-membership ðV new Þ of X t by using centres of clusters (v Ã Þ obtained in Step 1.
Step 7. Apply lag process for U new and V new obtained in Step 6 and obtain new input matrix ðX new input Þ; the new standardized input matrix ðZ new input Þ and standardized input matrix of the test set (Z test input ) that includes z scores of lagged variables of memberships and non-memberships.
Firstly, the new input matrix ðX new input ) that includes lagged variables of memberships and non-memberships is obtained by using U new and V new : Then, elements of X new input matrix is transformed to standard normal values, and so, Z new input matrix has obtained that means and standard deviations of each column of Z new input are 0 and 1. Lastly, Z test input that includes the last some rows of Z new input is obtained. For example, X new input and matrix obtained by using U new and V new obtained for X t ¼ f8; 10; 11; 12; 16; 13; 14g, c ¼ 3 and m ¼ 2 is given below: In Eq. (35), Z test input is the last row of Z new input given in Eq. (34) because n test ¼ 1.
Step 8 Step 4. Also, important principal components are selected from principal components by taking into consideration the criterion that the explained variance ratio needs to be at least 0.85. In this way, IPC test matrix is obtained. Lastly, the new input matrix ðP test input Þ is obtained by putting a vector that consists of 1 value to the first column of IPC. Step 9. Calculate the forecasts of X test t ð b X test t Þ by using b b and P test input obtained Step 5 and Step 8 via Eq.
For example, b X test t calculated by using b b and P test input obtained in Eqs. (30) and (39) are given below for Step 10. The root of mean square error (RMSE) and mean absolute percentage error (MAPE) values are calculated by using b X test t obtained in Step 9 via Eqs. (42) and (43). Some of the distributions of the 25 PEC time series used in the application are similar and the others are different. When the 24-time series graphs are viewed in general, it can be seen that the 24-time series may have characteristics Fig. 2 Yearly PEC observations of some countries in the million tonnes of oil equivalent such as linear increasing trend, curvilinear increasing trend, and stationary or random fluctuations. Some graphs of 24 PEC time series that differ from each other are given in Fig. 2.
The proposed method is compared with the following benchmarks.
DN-PSO Method: The Dendritic neuron (DN) method is an artificial intelligence model proposed by Zhou et al. (2016)  ARIMA Method: Autoregressive integrated moving average (ARIMA) method is a classic time series model that is often preferred by many researchers for time series analysis Kadılar (2005). Box-Jenkins's (1976) approach is an effective method that provides to be obtained the best ARIMA model. A Matlab m-file'' is prepared for solving PEC time series with Box-Jenkins's (1976) approach.
H-ES Method: Holt's exponential smoothing (H-ES) method is one of the exponential smoothing methods for time series analysis Kadılar (2005). A Matlab m-file'' is prepared for solving PEC time series with the H-ES method.
A-ES Method: Ata's exponential smoothing (A-ES) method proposed by Yapar (2016) is a new exponential smoothing method in the time series literature. A Matlab m-file'' is prepared for solving PEC time series with the A-ES method.
How applications are made with the proposed method is given below.
• Each of the 24 PEC time series has 52 observations.
The first 47 and the last five observations of the time series are used for the training set and the test set, respectively. In other words, the length of the test set is determined as 5 for all applications. • The number of fuzzy sets (c) is tried between 3 and 10 increasing 1 unit. • The degree of the model (m) is tried between 2 and 5 increasing 1 unit. • The architecture that has minimum RMSE value is selected as the best architecture, and so, the best forecasts are obtained for the test set. • 30 runs are executed using random values (rn), and so, some statistics are calculated for RMSE and MAPE values.
How applications are made with alternative method is given below.  Bas et al. (2015) were applied like those applications of the proposed method above-mentioned. • In applications of the SC-FTS proposed by Song and Chissom (1993a) method and C-FTS method proposed by Chen (2002), the number of fuzzy sets is tried between 3 and 10 increasing 1 unit. • The methods of the ARIMA method Box-Jenkins's (1976), Holt's exponential smoothing (Kadılar (2005)) and ATA method proposed by Yapar (2016) method were applied as classic time series methods.
In the applications, DN-PSO method proposed by Zhou et al. (2016), ANN-PSO method proposed by Yolcu et al. (2013), Median ANN method proposed by Aladag et al. (2014), Median-Pi-ANN method, FTS-N method and the proposed method were the methods that can be solved by trying different random initial parameters. For this reason, for each PEC time series, these methods were tried 30 times by using 30 different random seed values, and so, 30 different RMSE and MAPE values were obtained for each time series in the applications of this study. The mean and median of these RMSE values are given in Table 1. When  Table 1 is examined, it is seen that the RMSE statistics of the proposed method are smaller than the RMSE statistics of the other methods for the vast majority of 25 PEC time series solved. Therefore, it can be said that the proposed method is largely better than other methods in forecasting the future.
In the applications, for each of 25 PEC time series, minimum RMSE values of all methods used in this study are given in Table 2. When Table 2 is examined, it is seen that minimum RMSE values of the proposed method are smaller than minimum RMSE values of the other methods for the vast majority of 25 PEC time series solved. Therefore, it can be said that the proposed method has a largely better forecasting performance than other methods.
In addition, success rates obtained by examining Tables 1 and 2 were given in Table 3. When Table 3 is examined, it is seen that the proposed method has the highest success rates with 56%, 64% and 44% for minimum RMSE means minimum RMSE medians and minimum RMSE values, respectively. For instance, according to Table 3, 16 RMSE medians (16/25 = 64%) of 25 RMSE medians of the proposed method are smaller than RMSE medians of the other methods. Other methods have quite lower success rates than the success rates of the proposed method. Therefore, it can be said that the proposed method has a largely better forecasting performance than other methods according to RMSE measurement. Similarly, the proposed method has also the highest success rates with 44%, 44% and 32% for minimum MAPE means, minimum MAPE medians and minimum MAPE values, respectively. So, it can be said that the proposed method has a largely better forecasting performance than other methods according to MAPE measurement. However, the success rates with scores of 36% and 32% of ANN-PSO method proposed by Yolcu et al. (2013) for MAPE means and median means, respectively, and the success rate 24% of FTS-N method proposed by Bas et al. (2015) for minimum MAPE values are remarkable high values. Thereby, as a result, it can be said that the proposed method is the best, ANN-PSO method Yolcu et al. (2013) is the second-best

Explainable results of the proposed method
The inputs in the model can be tested and commented on. The inputs of the models are scores that are obtained from the principal component analysis. The contribution of each score to the forecasting model can be tested and commented on. An example for the Turkey time series is given below by using tables and graphs. In the forecasting model of the Turkey time series, six PCA score variables are employed as inputs. It is seen that all inputs are significant in the model in Table 4. The importance of the inputs can be determined according to coefficients. The importance graph is given in Fig. 3.

Conclusions and discussion
The high-order definition of intuitionistic fuzzy time series definition of Egrioglu et al. (2019b) means that any time series is affected by lagged time series, lagged membership values and non-membership values. However, indeed, membership and non-membership values are obtained from the lagged time series. So, it can be thought that membership and non-membership values include, largely, the information of lagged time series. From this viewpoint, in this study, a new high-order fuzzy intuitionistic time series definition given Definition 1 is introduced. It is explained that any time series is affected by lagged membership values and non-membership values in Definition 1. For this reason, it can be said that the advantage of the proposed method that uses Definition 1 is simplifying the model, since the number of inputs in the model is reduced. In the literature, intuitionistic fuzzy time series models have been proposed in recent years, since the definitions of classical fuzzy time series in literature can be insufficient for many real-life series. The reason for this is fuzzy time series approaches use only membership values for forecasting any time series. However, both membership and non-membership values are used in the intuitionistic fuzzy time series methods. Thus, the intuitionistic fuzzy time series methods are the approaches to increase the forecasting performance because of being used of more information. Also in this study, similarly, it has been shown that the forecasting performance of the proposed method that has a new intuitionistic fuzzy time series approach has much better forecasting performance than the classical time series and the classical fuzzy time series methods. Therefore, in future studies, it could be suggested that researchers prefer to use the intuitionistic fuzzy time series approaches instead of the classical time series and the classical fuzzy time series approaches.
In this study, a new forecasting algorithm is proposed, and it is shown that the proposed method has better forecasting performance with the best success rates between 32 and 64% than many time series models in the literature. In the proposed method, the principal component analysis method was used to make the dimension reduction. Then, the robust regression is used in the stage of determination relationships of the proposed method. Therefore, it can be suggested that some regression methods or statistical methods are used in the stage of determination relationships of the intuitionistic fuzzy time series that will be made in the future. In this study, the forecasting performances of 11 different time series forecasting methods are compared with each other. Then, for each PEC time series in the Europe-Eurasia continent, the best methods and the parameters of the best methods are given in Table 2. So, researchers that have worked in the energy fields could decide which methods they will select by looking at Table 2, when they would like to solve any PEC time series in the Europe-Eurasia continent.
In fuzzy time series methods, inputs of the models cannot be tested or commented on. The proposed method is an explainable artificial intelligence method. The proposed  Importance of Inputs Fig. 3 The importance graph of model inputs method is the first explainable fuzzy time series method in the literature. It is possible to test, order and comment on all inputs in the model because the fuzzy relations are modeled by using robust regression analysis. As a result of this study, it can be said that explainable artificial intelligence methods can be produced by combining statistical and machine learning methods. It can be expected that the studies about explainable fuzzy time series methods will be increased soon.

Declarations
Conflict of interest The author declares that there is no conflict of interest toward the publication of this manuscript.