A Collaborative Multiplicative Holt-Winters Forecasting Approach with Dynamic Fuzzy-Level Component

The adoption of forecasting approaches such as the multiplicative Holt-Winters (MHW) model is preferred in business, especially for the prediction of future events having seasonal and other causal variations. However, in the MHW model the initial values of the time-series parameters and smoothing constants are incorporated by a recursion process to estimate and update the level (LT), growth rate (bT) and seasonal component (SNT). The current practice of integrating and/or determining the initial value of LT is a stationary process, as it restricts the scope of adjustment with the progression of time and, thereby, the forecasting accuracy is compromised, while the periodic updating of LT is avoided, presumably due to the computational complexity. To overcome this obstacle, a fuzzy logic-based prediction model is developed to evaluate LT dynamically and to embed its value into the conventional MHW approach. The developed model is implemented in the MATLAB Fuzzy Logic Toolbox along with an optimal smoothing constant-seeking program. The new model, proposed as a collaborative approach, is tested with real-life data gathered from a local manufacturer and also for two industrial cases extracted from literature. In all cases, a significant improvement in forecasting accuracy is achieved.


Introduction
In the contemporary world of business, due to frequent technological advancements and operational changes, competition, complexities and challenges have increased drastically [1]. Moreover, easier access to international markets, expeditious technological advancement and the ceaseless requirements of quality products has complicated the scenarios for manufacturers. Hence, to contend with this dynamic but uncertain environment, manufacturers have no choice other than hastening their steps to gratify customers' demands on time and efficiently [2]. In this regard, the urge for a more accurate and reliable forecasting approach cannot be overlooked.
Plenty of forecasting approaches with innovative ideas have been proposed in the recent past which are broadly classified as qualitative and quantitative techniques. In the qualitative approach, prediction is made on the basis of the market dissection and the judgement of the experienced managers, whereas, the quantitative techniques mainly focus on historical data. However, in the latter forecasting class, the decision-making process is generally based on more or less complex mathematical models [3].
Among the quantitative approaches, researchers have shown particular interest in exponential smoothing techniques characterized by simplicity, robustness and ease of use [4]. Consequently, a number of exponential smoothing techniques have already been proposed, modified and/or improved in literature to anticipate an unknown event, enhance forecasting accuracy, as well as optimize a business plan e.g., production plan, production schedule and inventory control etc. [4][5][6][7][8][9]. Nevertheless, with this persisting evolution of techniques, the single and/or simple exponential smoothing model of Brown [10], double and/or trend correction model of Holt [11], triple and/or Holt-Winters approach of Winter [5] are still considered as reliable as well as the widely adapted forecasting approaches in literature.
Interestingly, all of these three exponential smoothing techniques are observed to follow the same principle in forecasting i.e., they estimate future events through the extrapolation of historical data patterns. However, with the presence of similarities in principles, these exponential smoothing techniques also demonstrate a number of distinguishing characteristics in their computational processes, parameter counts and responsiveness to different time-series components [9][10][11][12]. For instance, the simple exponential smoothing technique fits well with the time series exhibiting the changing mean having no seasonal and trend components. Whereas, with the presence of a trend component, varying level and growth rate, the Holt's trend correction model shows its superiority over the simple one. Meanwhile, the Holt-Winters forecasting model is observed to outperform the two former techniques for the time series, having changing seasonality, mean and growth rate [12][13][14][15].
Nevertheless, by featuring the mode of incorporation of these parameters i.e., either in a constant rate or by an incremental amount, the Holt-Winters model is further classified into additive and multiplicative models [9,[12][13][14][15][16]. Among these, the multiplicative Holt-Winters model, which was originally proposed by Winter in 1960 [5], is adopted for the time series in which the seasonal variations fluctuate along with the mean values [17]; whereas the additive model is chosen for the time series in which seasonal variations do not fluctuate along with the mean values [16].
However, to incorporate the effects of changing level, growth rate and seasonal pattern in forecasting, both of these Holt-Winters approaches practice recursive equations, the appropriateness of which depends largely on the initial values of the level, growth rate and seasonal components [14,18]. In this regard, although the initial values of the level and growth rate are usually estimated by using the least square regression method, the consideration of stationary mean adversely impacts the forecasting accuracy. However, it is computationally difficult and poses an even greater challenge to dynamically estimate and periodically update the values of the level (L T ) by modifying the conventional multiplicative Holt-Winters (CMHW) approach.
Hence, in this paper a collaborative forecasting approach is proposed in which the component level (L T ) is estimated by a developed fuzzy-based model and embedded into the CMHW approach. To do so, at the initial stage a simulation study is conducted to determine optimal values of the smoothing constants; whereas, in the second stage, the parametric values of the multiplicative Holt-Winters approach are determined and the parameter level (L Tpred ) is predicted through the incorporation of a developed fuzzy-based prediction model. The estimated level (L Tpred ) is then embedded into the conventional Holt-Winters method for each incremental time step to forecast the demand. The adequacy of our proposed collaborative forecasting approach as well as the level prediction model is assessed by the relevant metrics such as goodness of fit, mean absolute deviation (MAD) and mean squared error (MSE). Subsequent sections of this paper are arranged as follows: Section 2 presents a literature review on different demand patterns, time-series characteristics, forecasting approaches and the applications of artificial intelligence (AI) in forecasting. Section 3 highlights the research method, while Section 4 discusses the developmental stages of the fuzzy logic-based level prediction model for seasonal demand forecasting. The application of our proposed approach is described in Section 5. The corresponding results and discussions are provided in Section 6. Finally, in Section 7 relevant concluding remarks are made.

Literature Review
Efficient management of operations in manufacturing industries and services is important [19][20][21][22][23] and requires an accurate prediction of customer demand as this directly influences strategic, tactical and operational decisions. More specifically, a higher level of forecasting accuracy aids the manufacturer in formulating the appropriate production plan, restoring balance, scheduling the production line and maintaining proper inventory level as well as satisfying customer demand on time [24].
In the domain of forecasting, the time series is a set of values for one or more variables presented through consecutive time increments [25] and used in engineering, science, economics, finance, meteorology, telecommunication, etc. The hours, days and months are contemplated as the periods of the historical observation and/or the time series [26]. The forecasting methods featuring these time series are mainly categorized as regression, smoothing and decomposition techniques [27], among which, the exponential smoothing approaches are found to be the most widely used due to its robustness [9,28], simplicity [9] and the small number of required variables [2,9]. These favourable characteristics justifying the ease of implementation and cost effectiveness of this class of forecasting method [2,29,30] have attracted numerous researchers to dig continuously for further improvement [31].
Among these exponential smoothing approaches, the three most commonly adopted and addressed in literature can be named as simple exponential smoothing [10]; double exponential smoothing and/or Holt's method [11]; and triple exponential smoothing and/or the Holt-Winters method [5]. Simple exponential smoothing is renowned for its simplicity and ease of use, whereas Holt's method gives better results for the time series having trend components. Despite having superior performance for such cases, Holt's method loses its exactitude while dealing with seasonal time series [32]. However, in the presence of both the seasonality and the trend in time series, the Holt-Winters method shows better performance with minimal computational effort [24].
Two forms of Holt-Winters method (additive and multiplicative) are found to be the most widely adopted approaches in literature for seasonal demand forecasting, among which the additive forecasting method is observed to be suitable apparently for a situation where the seasonal peak adopts a consistent pattern throughout the demand horizon, and the latter fits well for a dynamic increase or decrease of the seasonal peak. However, during the implementation of these techniques, difficulties are encountered in determining the initial values of the time-series parameters (i.e., level, growth rate and seasonal index), quantifying the optimal values of smoothing constants, as well as using the intricate recursive computational processes [24,31,33]. As a consequence, and similar to the other time-series forecasting techniques, several approaches and modifications have also been proposed to improve the accuracy of the conventional Holt-Winters method by featuring the relevant issues.
For instance, Grubb and Mason in 2001 [34] proposed a modified Holt-Winters model for forecasting UK air passengers, where the values of the smoothing constants were determined by minimizing the mean square error (MSE). In the research of Lepojević and Anđelković-Pešić in 2011 [35], the MSE is also used as a performance measure and is a prerequisite for determining the appropriate trio of the smoothing constants. On the other hand, to estimate the initial values of the level, growth rate and seasonality, Chatfield and Yar in 1991 [17] considered these three parameters as the autonomous variables, and determined their initial values by minimizing one-step-ahead forecasting errors. Bermudez et al. in 2006 [24] also computed the initial values of these three time-series parameters by minimizing the forecast error. Meanwhile, Co and Boosarawongse in 2007 [4] adapted the Microsoft Excel solver tool for solving the mathematical interpretation of the Holt-Winters method where the initial values were determined by minimizing the forecasting errors. An analytical approach is also proposed by Bermudez, Segura and Vercher [31] to quantify and optimize the initial values of the time-series components. Moreover, Tratar in 2010 [36] also used the forecast error as an objective function and solved the mathematical interpretation of the Holt-Winters approach (through MS solver) by optimizing the initial values.
In addition, while forecasting the short-term electricity demand through the Holt-Winters approach, Taylor [37] considers the averages of the early observations as the initial smoothed values of the components. Arora and Taylor [38] also take the mean of historical data as the initial values of time-series components. Sudheer and Suseelatha [39] consider the average of the first day's historical data of their test week as an initial level and the average of the first two days' historical data as the initial growth rate. The initial seasonal components are calculated by averaging the yearly historical data. Meanwhile, Costantino et al. [40] average the historical data to estimate the initial values of the time-series components.
In order to forecast, the multiplicative HW method needs to undergo an intricate recursive procedure for updating its level, growth rate and seasonal index. The exactitude or appropriateness of this recursion process largely depends on the initial values and the patterns of the historical data [18,27]. Nevertheless, with the availability of various estimation processes for the initial values as presented earlier, the recent research works consider the stationary means in initiating the recursion process. However, incorporation of initial values in the stationary form is anticipated to affect forecasting accuracy in an adverse way. On the other hand, the computational difficulty of dynamically quantifying and updating the values of level (L T ) by modifying the conventional multiplicative Holt-Winters (CMHW) approach poses a greater challenge.
In these situations, an intuitional intervention within a mathematical technique known as a collaborative approach might lead to better performance [3,41,42]. In this context, fuzzy logic-based prediction models could aid in quantifying the levels and improving the overall forecasting performance. Nowadays, the fuzzy logic based model is widely applied in various fields of engineering and applied science [43][44][45][46][47][48][49]. For instance, to enhance the wastewater treatment process effectiveness by controlling the volume of dissolved oxygen in the reactor, Yang et al. [50] proposed an archetype of a fuzzy logic-based controller. This proposed predictive controller exhibits a satisfactory performance, usage of the conventional sludge model, and least square technique to define the rules and the parameters, respectively. Besides, the researchers have also found a way to employ the fuzzy logic in their proposed model for enhancing system stability and performance by means of fault detection [51,52].
Separately, by considering the need of developing an appropriate schedule for satisfying a customer's electricity demand, Mamlook, Badran and Abdulhadi [53] proposed a fuzzy inference system (FIS)-based approach to predict the short-term electricity load. In developing their intended fuzzy logic-based model, the authors have followed four successive steps i.e., fuzzification, development of fuzzy knowledge base, construction of the fuzzy inference engine and de-fuzzification of the outputs. Besides, the researchers have also considered five different input parameters, their corresponding membership functions and 37 different fuzzy IF-THEN rules in developing the Mamdani inference system. Thus, they were successful in reducing the mean absolute percent error (MAPE) by 5% in comparison to the traditional statistical forecasting approaches.
Cheikhrouhou et al. [41] proposed a fuzzy event-based collaborative forecasting technique by featuring the transient factors, transferred impact factors, quantum jump factors and trend change factors. In their proposed collaborative technique, the Auto Regressive Integrated Moving Average (ARIMA) model has been adopted as the initial mathematical approach for the inclusion of judgmental consideration. However, it is to be noted that the results obtained by their proposed collaborative technique outperforms those yielded by the ARIMA model. Meanwhile, Huang, et al. [54] also proposed an FIS for predicting the real-time electricity loads of smart homes. In their proposed approach, with four general steps of developing the FIS, the complexities of defining the associated locations of the membership functions are solved through a particle swarm optimization (PSO) algorithm. By using the developed FIS, the authors managed to reduce the MAPE compared to a third-order auto-regressive model and an artificial neural network approach.
However, despite the wider acceptance of the Holt-Winters and fuzzy logic-based approaches in forecasting, no attempt has been made yet to combine them in order to avoid the effect of stationary means on the forecasting accuracy. Hence, in this research endeavor an effort is made to propose a collaborative forecasting approach where the component level is estimated by a developed fuzzy-based level prediction model embedded into the multiplicative Holt-Winters approach for each incremental time step to safeguard higher forecasting accuracy.

Proposed Methodology
The approach presented in this research is intended to propose a collaborative Holt-Winters model for probable exclusion of the effect of stationary means and to make forecasting for seasonal demand more accurate. In this context, the values of appropriate smoothing constants are attained by simulating the historical demand data through an encoded program in MATLAB (R2009b, MathWorks, Natick, MA, USA) to be termed the smoothing constant simulator. Meanwhile, the initial values of the level and growth rate are estimated by the least square method. The initial values of the seasonal components are estimated by de-trending the time-series data, calculating the mean of the de-trended data for each period of chosen seasonal length, and then normalizing the obtained mean values to unity by using the correction factor [27].
All these values are used in determining the time-series parameters such as the level (L T ), growth rate (b T ) and seasonality (SN T ) as well as the forecast (Y T ) by the conventional multiplicative Holt-Winters method. Moreover, the determined values of three parameters would also help to estimate the error (E L ) that is perceived to be magnified by the stationary nature of the initial values. In the fuzzy logic-based level prediction model the parameters L T , b T , SN T and the actual demand Y T are taken as input variables, while the values of the predicted level (L Tpred ) is considered as the output. A Mamdani fuzzy inference system is applied to contemplate the judgmental knowledge and the relationships established between the variables by means of IF-THEN rules. A de-fuzzification method is employed to obtain the non-fuzzy levels (L Tpred ) which are embedded into the multiplicative Holt-Winters method to obtain the collaborative multiplicative Holt-Winters (Coll-MHW) model-based forecasted demand (Y Tc ). Finally, the accuracy of our proposed collaborative model is assessed by comparing the MAD, MSE and the coefficient of determination for real-life industrial data and two case studies extracted from literature. A brief summary of the research methodology is depicted in Figure 1 by a flow diagram.

Development Stages of the Proposed Approach
As stated earlier, the proposed collaborative forecasting approach consists of a number of successive steps such as introduction of the conventional multiplicative Holt-Winters method, definition of the optimal smoothing constants, determination of the parametric values and forecasted demand, development of the fuzzy logic-based model, and adjustment of the predicted parameter level. All of these steps are discussed comprehensively in the following subsections.

Development Stages of the Proposed Approach
The multiplicative Holt-Winters model for seasonal demand forecasting consists of an iterative process by incorporating the estimated level (L T ) defined as: Here, D T is the actual demand for period T, SN T−l represents the seasonal index for period T − l, where l denotes the seasonal length, and b T is the growth rate expressed as: In the Holt-Winters method, the seasonal factors are updated at each iteration. Thus, SN T and SN T−l , the seasonal indices respectively for the periods T and T − l, are updated by using the following equation: Generally, starting with the data of a specific period, the Holt-Winters method is able to forecast for the next one or more periods denoted as τ steps. Thus, for example, τ = 3 means that it is possible to forecast for three periods, from period T to T + 1, T + 2 and T + 3. The τ step ahead of forecasted demand (Y T+τ ) can be exposed as: So, one step ahead forecasted demand reflecting τ = 1 can be derived from Equation (4) and presented as Equation (5). The error, E T(L) can be determined by Equation (6). This equation represents a computational process to measure the deficit in estimating the level parameter for each incremental time step, T. In other words, the error, E T(L) expresses the de-seasonalized and de-trended one-step-ahead likelihood deficit which is used for training the rule base of the fuzzy inference system.
As mentioned earlier, the three smoothing constants (α, β, and γ), also known as the trios and used in the conventional multiplicative Holt-Winters method, are often ascertained subjectively or by using some iterative techniques. In this work, the optimal values of the trios are determined by simulation through a coded program in MATLAB.

Development Process of Smoothing Constant Simulator
To determine the optimal values of the trios (α, β, γ) as the smoothing constants within a range of [0, 1], having an incremental step of 0.1, an algorithm is formulated and a program coded in MATLAB, the main steps of which are stated below: The Pseudo code of Smoothing Constant Simulator is given as follows: Initialize parameters 3.
Calculate first step ahead forecast: 7.
END FOR 10.
END FOR 11.

Development Process of Smoothing Constant Simulator
In fuzzy logic-based prediction model the knowledge base, as shown in Figure 2, comprises two components-the data base and the rule base. Usually, the data base is represented by the formulated membership functions, whereas the fuzzy IF-THEN rules form the rule base of the intended knowledge base design. However, the core of this kind of knowledge-based model is the inference engine that employs a set of fuzzy IF-THEN rules to conduct the reasoning mechanism. Generally, the fuzzy reasoning mechanism is defined as a structured procedure for quantifying the conclusion from a given set of rules and conditions. Nevertheless, a number of fuzzy inference-engine modelling approaches are readily available in literature e.g., the Takagi-Sugeno-Kang (TSK), the Mamdani, etc. Among these, the Mamdani has advantages over others in terms of choosing particular fuzzification and defuzzification procedures, relaxation in representing the membership functions, as well as the wide range of logical operators [55]. Hence, by featuring computational and modelling flexibilities, the Mamdani system is adopted for our research work.
In the process of embedding the defuzzified values of level in the MHW model, computation begins with conversion of the numerical inputs into 'linguistic variables' and fixation of their ranges. Upon defining the linguistic variables for fuzzification, the membership functions are also formulated and allocated to each of the variables. By depending on the nature of the models different types of membership functions e.g., triangular, trapezoidal, Gaussian etc. are used in practice. The triangular membership function is adopted for our model due to its simplified features. However, for validation purposes of the developed model in the fuzzy logic tool box, Equation (7) is used to articulate the fuzzification process by representing the input-output membership functions in mathematical form [42,56].
In Equation (7), x represents the value of the input or output variables, whereas c 1 , c 2 , c 3 are the coefficients of the membership functions. Besides, to build the relationships among input-output variables, several, 'If-Then' logical rules are developed for the inference engine. Meanwhile, the fuzzy outputs are derived by aggregating the fuzzy inputs through the logical operators "AND", where the inference engine brings the compositional effects by incorporating the 'If-Then' rules. Then the non-fuzzy values of the desired output are obtained by Equation (8) [56,57].
where a i represents the position of the respective singleton in the respective universe and µ i stands for the firing strength of the truth degree of the respected rules.

Estimation of the Parameter Level (LTpred) and the Forecast
Similar to the procedure followed by Hossain et al. [57], the level (L Tpred ) obtained from the fuzzy logic tool box is validated initially by comparing with the values found mathematically. The validated values are then embedded in the conventional MHW method to compute the collaborative multiplicative Holt-Winters (Coll-MHW) model-based forecasted demand (Y Tc ) by applying Equation (9).

Accuracy Assessment
To judge the accuracy of the forecasting approaches, several benchmarking indicators such as mean absolute deviation (MAD) or mean absolute error (MAE), mean squared error (MSE), mean absolute percent error (MAPE), root mean square error (RMSE) are quite common and adopted in [58][59][60][61]. However, for assessing the strength of our proposed approach, the MAD, the MSE and the coefficient of determination were determined and compared with those obtained through the conventional Holt-Winters method. The equations for the three metrics are as follows [57,58,60,62,63]:

Application of Our Proposed Approach
As mentioned earlier, three case studies are adopted to demonstrate and validate the applicability of our proposed approach. These cases are based on (a) demand data for transformer tanks collected from a local manufacturer; (b) average sales data of a Slovenian company (i.e., Danfoss District Heating) as extracted from literature [36]; and (c) demand data for plastic bags of a manufacturer in south of Spain as extracted from the literature [41]. The whole process is discussed step by step in the following sections.

Case Study-1: Demand Data for Transformer Tanks
The demand data collected from a transformer tank manufacturing company are adopted to validate the proposed methodology. The company mainly produces three types of distribution transformer tanks such as 1000 KV, 750 KV and 500 KV. Among the three models, the 750 KV transformer tank is found to be the vital item for the company. Since the management reviews its production plan at the end of each month, the monthly demand for 750 KV transformer tanks is considered in this study. As presented in Figure 3, the historical records of monthly demand for transformer tanks over the past four years ranging from 2010 to 2013 were revealed by the company. The first three years' data are used for the training and/or prediction model development purpose while the remaining data are used to test the proposed approach. However, prior to forecast, it is common to examine a time series for characterizing the presence of different demand patterns such as the randomness, trend, seasonality and cyclical behavior. In this context, a scatter diagram of time series representing the demand for 750 KV transformer tanks as shown in Figure 3 reveals an overall uprising trend with multiplicative seasonality.

Determination of the Smoothing Constants
MATLAB coding is used to simulate and determine the best trio of smoothing constants. As a primary objective of the simulation, the best combination of trio was selected against the lower MAD values of the forecast. Optimized values of the trio (α, β and γ) are secured by converging the MAD to its lowest possible value. The converging behavior of MAD is demonstrated in Figure 4a Figure 4a presents the results of the simulation in which the MAD values are changed in a line with the variation of β and γ while α remains fixed at its optimal value (α = 1). Figure 4b illustrates the effect of α and γ on MAD values with a fixed value of the variable β (β = 0), whereas Figure 4c shows the converging behavior of MAD with respect to α and β for a constant value of 0.2 for γ. It can also be seen that in all of these three exhibits, the MAD is converging to its unique minimal value of 0.89. Hence, the optimum values of the trio (α, β and γ) as obtained through simulation are 1.0, 0.0 and 0.2, respectively.

Determination of the Initial Values
Once the optimum values of the trio (α, β and γ) are obtained, the level (L T ), the growth rate (b T ), the seasonal parameter (SN T ), the forecasted demand (Y T ), and the error (E T(L) ) are computed respectively by using Equations (1)- (6). The computed values are presented in Table 1.

Development of Level Prediction Model
In practice, the Holt-Winters method has been recognized to have relatively higher accuracy and reliability compared to other conventional forecasting approaches. However, the potential lack caused by the consideration of stationary initial values in its recursion process has left scope for incorporating artificial intelligence to improve the forecasting accuracy. In this context, a level prediction model is built via a fuzzy logic-based model. Salient steps in the procedure are described in the following sections.

Choice of Inference Mechanism and Variables
In this research, the developed fuzzy logic-based model is implemented in the Fuzzy Logic Toolbox of MATLAB R2009b where the Mamdani inference system is chosen for predicting the parameter level (L Tpred ). It is also necessary to define the input and output variables since their interrelations will materialize the objective(s) of the fuzzy logic-based model. Input variables i.e., the demand (D T ), the level (L T ), the growth rate (b T ), the seasonal parameter (SN T ), and output variables i.e., the predicted level (L Tpred ), in the MATLAB platform are exhibited by Figure A1 in Appendix A.

Fuzzification and Formulation of Membership Functions
In the fuzzification stage, the values of all the linguistic variables, as shown in Table 2, are fuzzified with the aid of membership functions. Usually the membership functions are considered as the core of fuzzy logic-based modelling as they have a direct influence on the strength and weakness of the model [64,65]. Among the different types of membership functions, the triangular one is selected for the fuzzification purpose due to its ease of construction and versatility for system representations [66]. However, in this intended development of the fuzzy logic model, the fuzzification of the linguistic parameters is conducted by the following functions: where, i 1 , i 2 , i 3 , i 4 , o 1 are the variables related to the level, growth rate, seasonality, actual demand and predicted level, respectively. The triangular fuzzy numbers for the input variables i.e., the demand (D T ), level (L T ), growth rate (b T ), seasonal parameter (SN T ), and the output variables i.e., the predicted level (L Tpred ) are given in Table 3 by characterizing the definition as well as the properties presented in the work of Chu and Tsao [67].   However, to demonstrate this fuzzification process through mathematical interpretation [42,57,68], the membership functions for the input and output variables are worked out by using Equation (7). Figure A2a-e as provided in Appendix A illustrate all the developed membership functions in the MATLAB platform. However, for fuzzification and/or articulation of the membership functions, the associated values of the coefficients of the entire input and output variables are shown in Table 4.  Meanwhile, as adopted by [69], the membership functions and the linguistic fuzzy sets of a particular input L T are also deduced by using Table 4, Equation (7), Equation (13) and presented through Equations (18)- (24). The membership functions for the remaining input and output variables are also developed in the same manner as shown by Equations (18)- (24).

Development of Fuzzy Inference Rules and Process
For Case-1, 30 decision rules, as presented in Figure 5, are formulated and inserted in the fuzzy logic tool box. Each rule has two distinct parts: premise and conclusion. The input variables are interpreted with the premise part as "if" whereas the output variables are taken with conclusion as "then". To clarify how a rule works, an example is given for one of the decision rules along with the interpretation: Rule 1: If L = E. low and b = V. low and SN = V. low and D = E. low, then L (pred) = L2. This means that if the level is E. low, growth rate is V. low, seasonality is V. low, demand is E. low; then the predicted level is L2. After formulating the decision rules and membership functions, the inference process is carried out. In the first step of the inference, the conclusion of each decision rule is derived for a particular input by a comprehensive comparison among the ranges of input variables. At this step of implementation, the logical operator 'AND' is used to connect the inputs. The overall conclusion is made by the aggregation of the selected decision rules to a single fuzzy set. This aggregation is carried out by selecting the truth degree of each rule through the min-max operation. Finally, the non-fuzzy values are obtained by de-fuzzifying the fuzzy sets through Equation (8).
To perceive the defuzzification process mathematically and assess the validation of the developed level prediction model on MATLAB [56,57], an example is presented with an input for L T = 16.74, b T = 0.27, SN T = 1.13 and D T = 20.00 corresponding to the 19th month. The membership functions, the decision rules as well as the rule viewers of the fuzzy logic tool box (as shown in Figure 6) indicate that for the given inputs only the decision rule 11 is fired (i.e., output for this rule is non zero) and all other rules remain in the non-fired (i.e., output of these non-fired rules are zero) position. Similar to what is adopted by Hossain et al. [69], for this particular rule the fuzzy values of the input L T , b T , SN T and D T are computed by using Equation (7) and the values obtained are as follows: Whereas, for rule 11, the firing strength (truth value) are determined as follows: = min(0.260, 0.644, 1.350, 1.00) = 0.260 Finally, for the defuzzification, the crisp output of the predicted level is computed as 17.34 by using Equation (8); whereas the developed fuzzy logic-based model in MATLAB as shown in Figure 6, provides an output of 17.4 for the same inputs. This shows that the obtained output of 17.34 is very close compared to the output 17.4 given by the model developed in the MATLAB fuzzy toolbox. By featuring this higher accuracy, the obtained output of the predicted level (L Tpred ) is used to compute the new series of forecasted demand (Y Tc ).

Prediction of Level and Consequent Forecast
At this stage, the new series of collaborative model-based forecasted demand (Y Tc ) for the test set is obtained by using Equation (9). The value of the demand (D T ), level (L T ), predicted level (L Tpred ), forecasted demand (Y T ) and the consequent model-based forecasts (Y Tc ) for the test period are given in Table 5.

Case Study-2: Average Sales Data of Danfoss District Heating
The average sales data of 4247 products of a Slovenian company named Danfoss District Heating, which was originally presented by Tratar [36], is adapted in this research for validating and assessing the strength of our proposed collaborative approach. The time series, as shown in Figure 7, represents the monthly average sales data for the time span of four years ranging from 2005 to 2008. The data set also depicts the varying means and the overall uprising trend with multiplicative seasonal patterns. However, to assess the strength of our proposed collaborative approach, initially the starting values of the level, the growth rate and the seasonal components are estimated. These initial values are incorporated with the conventional multiplicative Holt-Winters method to determine the time-series parameters such as the level (L T ), growth rate (b T ), seasonality (SN T ) as well as the forecast (Y T ). A fuzzy logic-based level prediction model is also developed to attain the output level (L Tpred ). To do so, the Mamdani fuzzy inference system is chosen to contemplate the judgmental knowledge and the relationships established between the variables by means of 39 IF-THEN rules and the logical operator 'AND' as shown in Figure 8. Finally, the de-fuzzification is performed by using Equation (8) to obtain the non-fuzzy levels (L Tpred ) which are embedded into the multiplicative Holt-Winters method to get the Coll-MHW model-based forecasted demand (Y Tc ) through Equation (9).

Case Study-3: Demand Data for Plastic Bags
In their research, Cheikhrouhou et al. [41] presented a four-year monthly demand data set for plastic bags collected from a manufacturer in the south of Spain. The data set as presented in Figure 9 exhibits the changing means, growth rates and seasonal patterns for which the multiplicative version of the Holt-Winters approach looks to be a relevant for forecasting. From this data set, the first three years' data (for 2004-2006) are used in the training period to forecast demand for the period of 2007. To forecast through our proposed collaborative approach, the starting values of the level, the growth rate, the seasonal components and the forecast are estimated initially by the conventional multiplicative Holt-Winters method. In addition, a Mamdani type fuzzy logic-based level prediction model is developed to attain the output level (L Tpred ) through incorporation of 47 IF-THEN rules and logical operator 'AND', as shown in Figure 10. Meanwhile, the values of levels (L Tpred ) and forecasted demand (Y Tc ) are obtained by using Equations (8) and (9), respectively. Finally, the accuracy of the forecasted demand for the period of 2007 is assessed by means of MAD and MSE.

Results and Discussions
The forecasted results obtained by means of the Coll-MHW are used to evaluate the MAD, MSE and coefficient of determination as the statistical measures for accuracy assessment.

Evaluation of the Forecasting Accuracy
To evaluate the forecasting accuracy, measures of performance such as MAD and MSE are used. In this context Equations (10) and (11) are used to determine the values of MAD and the MSE respectively for all the three cases. As presented in Table 6, the MAD and MSE obtained for the forecast values yielded by the proposed Coll-MHW method outperform all other methods. It is to be noted that MAD and MSE values for Case-1 are evaluated by using the conventional multiplicative Holt-Winters method (CMHW) whereas for Case-2 and Case-3 the values of MAD and MSE are obtained from literature. As stated in the published papers, the forecast values for Case-2 were determined by CMHW and improved multiplicative Holt-Winters method (IMHW), while those for Case-3 were evaluated by ARIMA and collaborative adjusted (Coll. Adj.) models.
It would also be interesting to demonstrate the forecasting performance of the proposed approach for each of the individual periods of time series. In this respect, Figure 11a-c portray a comparative scenario for the three cases. The comparative scenarios as presented in Table 5 and plots in Figure 11a-c consistently reveal the superior performance of our proposed collaborative approach in comparison to the conventional MHW, improved MHW, ARIMA and Coll. Adj. forecasting models.

Assessing the Coefficient of Determination
Apart from assessment of forecasting accuracy, the coefficient of determination is evaluated to reflect the strength and validation of the models. To check the appropriateness of our proposed model, the values for the coefficient of determination are computed and compared with those for other methods. For this purpose, the transformer tank demand over a period of one year, the average sales of Danfoss District Heating, and the data on plastic bag demand are used as before. Table 7 and Figure 12a-c exhibit the coefficient of determination. It is to be noted that the coefficient of determination computed by our proposed Coll-MHW approach is compared with those yielded by the CMHW for Case-1 and Case-2 and by the Coll.
Adj. for Case-3. Regrettably, it could not be assessed for the ARIMA and IMHW methods due to the absence of the required data. However, as delineated in Table 7 and in Figure 12a-c, the coefficient of determination (R 2 ) values validate the superior performance of our proposed Coll-MHW approach.

Conclusions
Irrespective of whether a firm is oriented to services or manufacturing, many operational decisions are highly dependent on the predicted values of certain events or issues. In this context, the role of an appropriate forecasting approach cannot be ignored. Although a prediction can never be perfect, an endeavor is always under way to search for an accurate forecasting approach. The developed collaborative multiplicative Holt-Winters (Coll-MHW) approach as proposed in this paper is facilitated by a fuzzy logic-based prediction model to dynamically estimate the L T value and embed it into the conventional MHW method. This newly proposed Coll-MHW approach is tested through several case studies and is found to be robust with superior performance. Moreover, our proposed approach does not require the inclusion of any judgmental consideration like transient factors, transferred impact factors, quantum jump factors and trend change factors as quoted by Cheikhrouhou et al. (2011) [41]. So, the proposed Coll-MHW approach is expected to be a potential forecasting method to address multiplicative seasonal demand forecasting problems and to secure its place as an effective tool for prediction. In future, the proposed method is likely to be enriched in an extended, robust form for wider application through the inclusion of lifecycle demand data and AI-like neural network. In this context, the impact of the appropriate choice and adoption of different kinds of membership functions on the performance metrics for this proposed collaborative forecasting model is also kept in mind for an extended scope of future research.