Road Freight Transport Forecasting: A Fuzzy Monte-Carlo Simulation-Based Model Selection Approach

: As important as the classical approaches such as Akaike ꞌ s AIC information and Bayesian BIC criterion in model-selection mechanism are, they have limitations. As an alternative, a novel modeling design encompasses a two-stage approach that integrates Fuzzy logic and Monte Carlo simulations (MCSs). In the first stage, an entire family of ARIMA model candidates with the corresponding information-based, residual-based, and statistical criteria is identified. In the second stage, the Mamdani fuzzy model (MFM) is used to uncover interrelationships hidden among previously obtained models criteria. To access the best forecasting model, the MCSs are also used for different settings of weights loaded on the fuzzy rules. The obtained model is developed to predict the road freight transport in Slovenia in the context of choosing the most appropriate electronic toll system. Results show that the mechanism works well when searching for the best model that provides a well-fit to the real data.


INTRODUCTION
Technological forecasting as a dominant force to create social change is a new gate to the future. This is because it arises from the interplay of public policy, strategic policy, and scientific methods to incorporate possible future technological development in management decision making. The ability to accurately predict the future is essential for decision processes in many scientific areas and real applications, let it be the planning, organizing, strategy formulation, policy-making, and so on [1,2]. Thus, particularly in the last few decades, time series (TS) forecasting has become a stimulating research matter with mounting attention to various areas [3][4][5][6] The capability to provide accurate forecasts is also significant in the field of broader planning of any essential road transport projects and investments. In this context, the planning of potential strategies must also take into consideration the accurate forecasting of future trends in freight transport. Such kind of forecasting is crucial in the case of developing any of future strategic policies on the state level, such as, for example, a design of a suitable tolling system on highways, where the freight traffic is usually prevailing, especially the transit one.
The study in this paper represents a continuation of the previous research of the authors of this paper (see [7]). There, the research was recognised as one of the essential studies conducted to select the most appropriate technology for the Electronic Tolling System (ETS) based on providing accurate forecasts of long-term future trends of Slovenian Highwaysꞌ freight transport. Nevertheless, contrary to the previous research, this paper is more methodological than the applicative one. Also, the main focus is more dedicated to a model selection than long-term forecasting. In the modelling process, the model selection step is one of the essential steps conducted during the model identification, which includes the determination of the modelꞌs structure and estimation of its parameters [8].
The main aim of the paper is to present a novel twostage fuzzy-based model selection procedure, which is combined with the Monte-Carlo simulations. The proposed mechanism represents an alternative to the commonly used model selection approaches, where the best model is chosen based on observing, e.g., the Akaikeꞌs AIC information criterion or Bayesian BIC criterion [9,10]. In some cases, these criteria do not identify the most parsimonious model, thus they have certain deficiencies that our approach strives to overcome.
In the first stage, for different combinations of orders p and q, the whole family of 1,..., i N  plausible ARIMA model candidates is generated with different levels of fit to the given freight transport historical data  , 1,..., For each of these N candidates, the J various information, statistical, and residual-based criteria verify this assumption, the 1 for each i-th ARIMA model candidate. From these errors, a further analysis clearly isolates the "optimal" index * i that belongs to the best ARIMA model. The designed FM-MCS model selection mechanism was first tested for the simulation example and then used in the real case of freight transport data. Prediction results show that the developed forecasting model can provide a well-fit to the real data.

LITERATURE REVIEW 2.1 Forecasting the Road Freight Transport
When it comes to the forecasting of traffic flow, the methods can be, in general, divided into two main groups. The first group includes short-term forecasting (STF) methods. They are mainly used for advanced traffic management in the sphere of intelligent transportation systems at the level of urban metropolitan areas. These methods can be further divided into the classical methods (e.g., the ARIMA and other Box-Jenkins models, Kalman filters, etc. [7]) and the artificial intelligence/machine learning/data mining methods (e.g., the neural networks, support vector machines, etc. [7]).
The second group includes the methods for long-term forecasting (LTF) and is mostly used for forecasting the goods' movements on the state level. Some of these methods detected in the literature are [7]: The Dynamic Factor models, ARIMA models, ARMA-GARCH models, neural or neuro-fuzzy models, support vector regression models, etc. A more in-depth review of the methods developed for the road freight transport forecasting can be investigated in the literature (e.g., see [7]).

Forecasting, Regression and Fuzzy Logic
The classical linear regression models often suffer from a modelling error and measurement noise in the addressed data, which induces a sort of fuzziness around the precise exact data [12]. Moreover, the linear regression models often cannot deal with different nonlinearities, complex uncertainties, and chaotic behaviour of observed stochastic processes effectively enough.
To be able to deal with highly nonlinear patterns in the observed data, it is often recommended to use hybrid models or to combine several different models [13][14][15]. Furthermore, to deal with the modelling error and measurement noise in the data that are "captured" in the fuzzy environment, Tanaka and his colleagues [16,17] have recommended a fuzzy regression approach in their pioneer works in the Eighties, where an interval prediction model has been conducted. Since the contributions [16,17] occurred, the great body of literature followed dealing with the fuzzy logic (FL) applied to the time series regression, analysis, and forecasting (see [12] for contributions in the previous century).
In this century, Tseng and his colleagues [12] have combined the ARIMA model and fuzzy regression model into the FARIMA model for forecasting the exchange rate between the NT dollars and US dollars. This model was based on the triangular fuzzification of crisp ARIMA parameters and has achieved noticeably better forecasts than the competitive classical ARIMA model. The extension of the Tsengꞌs work was done by [18], where the combining of the SARIMA (p, d, q) (P, D, Q) model with the fuzzy regression model was carried out resulting in the FSARIMA model. The model was used for forecasting the production value of Taiwanꞌs machinery industry and the sale of soft drinks. Last but not least, the contribution [19] 1 SCAN -Smallest canonical correlation method 2 ESACF -Extended sample autocorrelation function 3 CCA -Canonical correlation analysis 4 CVA -Canonical variate analysis has employed the Neural networks and fuzzy logic to model the nonlinear patterns of residuals of the ARIMA model to forecast various exchange rates and the gold price. Some other essential works from the field can be found in the literature (e.g., see [20][21][22]). To summarize, the fuzzy logic has been extensively used to improve the deficiencies of the classical forecasting models. On the other hand, there can be practically no study found that would involve the FL for the model selection purposes in the way as it was done in this paper.

Model Selection Criteria and Procedures
A lot of the single criteria for time series model selection have been developed besides the aforementioned AIC and BIC, such as, for example, the AICc (for small sample size), the Hannan and Quinnꞌs HQ, the minimum description length (MDL), the Kullbackꞌs information criterion (KIC) and its corrected version (KICc), to mention some of them (see [23,24] for details). However, when we are dealing with more complex TS models, such as ARMA models, it is in most cases relatively difficult to choose the appropriate order of AR and MA parts.
The monograph work by [23] reveals a wide-ranging review regarding the TS model identification and selection, including the ARMAX, ARIMAX, and other transfer function models. Here, the linear innovation regressionbased estimation methods, such as the HR method, and the KP method, are discussed. These methods can determine the model's order as well as its parameter estimates, while some of them can also be applied to the Transfer function models [23]. Choi in his work also introduces the pattern identification methods, such as the SCAN 1 method, the R and S array method, the GPAC method, the Corner method, the ESACF 2 method, and other. Another important class of model selection and identification methods are Subspace identification regression-based methods, which are specially designed for ARMAX modelsꞌ state-space representations. The main protagonists of these methods are the CCA 3 method and the CVA 4 method, proposed by [25], as well as the N4SID 5 method [26], and the MOESP 6 method [27].
Besides, Gomez and Maravall have proposed an algorithmic automatic procedure, which has been later also implemented in programs TRAMO and SEATS [28]. More recently, [29] proposed a state space framework for automatic model selection and forecasting based on exponential smoothing methods. On the other hand, authors [30] outlined the so-called RETINA 7 automatic model selection procedure consisting of four main stages (data building and sorting, isolation of candidate models, the search strategy, and a selection of the most appropriate model). There also exist several other model selection strategies and approaches, such as the strategies with Genetic and Evolutionary algorithms (e.g., [31,32]), or Branch and Bound strategies, to remark some of them. Contrariwise to all these approaches, our approach is quite different since it is grounded on the aforementioned FM-MCS based mechanism.   have the form: When dealing with the   ARIMA , , p d q models, the requirements of stationarity, stability, and invertibility must always be fulfilled according to the rigorous requirements of the Box-Jenkins modelling approach [8].
The invertibility implies that when a modelled time series is expressed in the form of the function of past observations, their weights should decline as one moves further into the past (i.e., the more recent observations should be taken into account more heavily than those from the more distant past).

Fuzzy Logic and Mamdani FIS Systems
Fuzzy logic offers tools to incorporate uncertainty and the human way of reasoning, thinking, and perception [11]. In fuzzy logic systems, the values of variables are articulated by linguistic terms such as for example, "large", "medium", "small", etc. Here, the relationships are expressed in a form of "if-then" rules (if -antecedent, else -consequent), while the fuzzified outputs and can be converted to "crisp" values by using the defuzzification methods [11].
When doing fuzzy modelling, the fuzzification is a process for defining the degree of the membership of a specific value being a particular fuzzy set. This can be done by evaluating the so-called membership functions (MF) belonging to the given fuzzy sets. The membership functions determine the degree of the membership of an individual element in a fuzzy set. Regarding the fuzzy rules, they are the basic elements for apprehending knowledge in fuzzy models [11].
The Mamdani and Sugeno types of fuzzy systems are the most common basic systems when applying the fuzzy logic model. The Mamdani system addresses the process states by using the linguistic variables as a means to control the fuzzy rules. It contains the following basic elements: A fuzzifier to convert the values of input variables to input fuzzy sets; the knowledge base; the inference engine, and the defuzzifier to transform the output fuzzy set to the crisp values of the output variable (c.f. Fig. 1). The main goal of the fuzzy inference engine is to make specific calculations with the fuzzy rules. Moreover, it is based on the following three foundations: aggregation, implication, and accumulation. Further details about fuzzy logic and Mamdani systems can be found in the literature (e.g., see [11,33]). Fuzzy set can be mathematically presented as [33]: is the membership function that determines a degree of the membership of elements x U  to the fuzzy set A. The MF ( ) A x  can be expressed as follows:

Figure 1 The fuzzy inference system
The value of the linguistic variable can be expressed by using the following triple: the name of the variable,   A X is the set of values it can take, while the "Rule" determines the relationship between the linguistic terms and their physical meaning. According to Sivandam et al. [33], the fuzzy rules are a group of linguistic statements that define how the fuzzy system should carry out a decision considering classifying the inputs or controlling the output. Moreover, in general, the fuzzy rules can always be written in the following way (for the Multiple-Input-single-Output (MISO) systems) [11,33]:  Fig. 2 shows the conceptual framework with a novel two-stage FM-MCS model selection procedure (see also Fig. 1 and Eq. (1), (2) and (5)). Since the measured freight transport data   ORIG , 1,..., y t t T  from 1990Q1 to 2019Q1 had contained some inconsistencies, errors, unusual spikes, abnormal and missing data, as well as outliers, the pre-processing and data cleaning was needed in order to improve the adequacy of the data (see Fig. 2, block B). This way, a consistent and error-free time series  , 1,...,

THE DESIGNED FM-MCS MODEL SELECTION MECHANISM
was obtained in tonne-kilometres (Tkm).
The TS was given in the cumulative and aggregated form obtained on the quarterly basis, which has given T = 117 time samples from the first quarter of 1990 to the first quarter of 2019 (see [7] for details).

Figure 2
The Conceptual Framework with a novel two-stage fuzzy-based model selection procedure The unit root testing was applied in the next step to examine the stationarity of the time series (c.f. Fig. 2, block B). For this purpose, two tests were accomplished, that is the KPSS (Kwiatkowski-Phillips-Schmidt-Shin) test and the ADF (Augmented Dickey-Fuller) test [8]. The tests had revealed that the ARIMA model is needed in the further modelling procedure. Also, it was discovered that the first order of differencing (i.e., d = 1) of the TS is satisfactory to achieve its stationarity. Thus, in the further procedure, one was faced with a differentiated time series

The First Stage of the Algorithm
The 1. stage in Fig. 2 The other 20 models did not pass this modelling step due to different issues and violations regarding the Box-Jenkins methodology (the none-invertibility, the models' error had not been a white noise, etc. [8]). Because of similar reasons, the higher orders of p and q were also not taken into consideration due to none-invertibility, illconditioned behaviour, and other similar issues. While processing block C, for each model candidate, besides the estimated parameters ˆi   andˆi   , the estimated forecasts   , y t i and ARIMA modelling errors were also computed (see Fig. 2, block F).
In the next step of the 1. Stage (see blocks D, E, G, K in Fig. 2), the J = 7 information, statistical, and residualbased criteria  

The Second Stage of the Algorithm
While designing the FM of the type MFIS (see Fig. 1),

FR s s S 
were established (c.f. Fig.   1, block G, Fig. 2, block J, and Eq. (5)), while simultaneously, the appropriate membership functions for all FM inputs and FM output were also constructed (c.f. Fig. 1, blocks F and H, and Eq. (5)). Regarding the fuzzy rules, to simultaneously satisfy demands 1 and 2, the main principle was that for low values of FM inputs, the FM output should be high and vice versa. The same logic was adapted for all the intermediate situations. Such kind of thinking has given several dozens of different fuzzy rules. Considering the MFs, for both the FM inputs and the FM output, the Gaussian membership functions were chosen with three possible classes: "low", "medium", and "high", [11]. When the fuzzy rules and MFs were constructed, the next step was to prepare the setting for weights , p q are the optimal orders of the AR and MA parts, while * i is the optimal index). This manoeuvre was based on assumption that the best ARIMA model is "hidden" behind that index  In other words, this can also be implicitly understood by the fact that the best ARIMA model at given weightsꞌ combinations most often has had a minimal FMs error compared to other suboptimal ARIMA models. A further mechanism for determining the optimal index * i of the best   * * * ,1, ARIMA p q i model can be investigated by carefully observing the Pseudo-code find_the_best_model() in Fig. 3 (see also Fig. 2). The pseudo-code is intentionally given in the un-optimised form for easier understanding.

Some Details of Testing of Designed FM-MCS Model Selection Mechanism for a Simulation Example
In order to verify the capabilities of the developed mechanism, it was first tested for the case of simulation example. For this purpose, the simulated  

PRACTICAL NUMERICAL RESULTS
The modelling process and all other calculations were carried out in the MATLAB technical computing environment, with the additional use of the three MATLABꞌs toolboxes, i.e., the Fuzzy Logic Toolbox, the Econometrics Toolbox, and the Statistics and Machine Learning toolbox.

The Results of the 1. Stage of Model Selection Process
When all steps of the 1. Stage (see blocks C, D, E, G, F, K in Fig. 2)

The MFs and Fuzzy Rules Obtained by the FM-MCS Mechanism in the 2. Stage
In the 2. Stage (see Fig. 2  j j j X X X j  , as well as for the FM output ( L M H , , Y Y Y for "low", "medium", and "high" class of the output MF for the FM output Y).A simplified idea behind the setting of the 1,..., s S  fuzzy rules can be illustrated in the following excerpt.
When observing the rules given in Eq. (7), one can see that for low values of the FM inputs (i.e., L is j j X X ), the FM output Y should be high (i.e., H is Y Y ), and vice versa. The same logic holds for all the intermediate situations. As it turned out, a total of S = 174 rules was enough big number that the just explained logic had been satisfactorily covered.

The Best ARIMA (4,1,2) Model Obtained by the FM-MCS Mechanism after the Completion of MCS Iterations (2. Stage)
The MCS procedure in Fig. 2 Accordingly, taking into account Eq. (1) and (2), the following form of the structure of the best model appears: ,ˆ1   time-samples (app. 76% of all data -the estimation/insample interval) were used to carry out the model selection procedure and to identify the best model. The other (last) 28 time-samples were used to test the predictive power of the model (app. 24% of the data -the test/out-of-sample interval). As can be seen from Fig. 8, the model provides surprisingly good predictions * ( , ) y t i of the real road freight transport time series ( ) y t (measured in Tones km). There can be noticed some very sophisticated details of the TS dynamics, which were not entirely captured by the model. The reason is that the freight transport data dynamics likely contains a complicated nature with perhaps even a certain kind of nonlinear behaviour included. Despite this, the model generally provides an encouragingly good fit to the real data, particularly considering the key movements of the TS trend.  The modelꞌs predictive performance should also be analytically verified by calculating the achieved measures of the quality of the derived forecasting model. Besides the well-known classical measures (e.g., the aforementioned MAE, MSE, AIC, BIC, etc.), there have been many new measures applied in the last decade, which had become significantly attractive during the well-known Makridakis competitions [34,35]. Tab. 2 depicts the meaning of some of these new applied measures (see [34,35] for corresponding formulas).

The Predictive Performance of the Best ARIMA Model
Measures such as (sMAPE, MAAPE, MASE, etc.) can be treated as the most reliable and were also used to check the modelꞌs quality. Tab. 3 shows the calculated measures for the case of the best ARIMA model in order to verify its quality and performance. Some of achieved values in Tab. 3 are compared with the prescribed thresholds recommended in the literature. Computed measures confirm the adequate quality and performance of the derived model.

Discussion and Practical Implications
Despite the fact that designed FM-MCS mechanism that gives the best ARIMA model works surprisingly well and provides a well-fit to the real road freight transport data, its prototype is currently developed for short to medium term forecasting only. Thus, similarly as in the previous research [7], it might have been appropriate to combine the FM-MCS algorithm with a Monte Carlo Scenario Playing (MCSP) mechanism for generating future scenarios about the road transport trends. The MCSP mechanism might have represented the basis for calculating the long-term interval-type forecasts, for which the corresponding probability intervals in future timepoints could be computed. Besides, there are likely many other practical examples and cases not only in the infrastructure and road transport planning sphere, but also in many other theoretical and practical fields, where the designed FM-MCS mechanism might have offered an opportunity to the scholars and practitioners to facilitate the model selection procedure.
In many real applications, as much accurate forecasting as possible is of paramount importance. Namely, even slightly inaccurate forecasts can lead to wrong strategic decisions and wrong financial investments, which can have catastrophic financial consequences, such as inappropriate investments in infrastructure, bad choice of expensive management and control systems, and the like. On the one hand, our innovative approach to modelling and choosing the best model is slightly complex, it is true. However, on the other hand, the latter is almost a fully automated generic mechanism that helps to extract the best forecasting model, which could be deployed with possible slight modifications for other similar real-world cases of forecasting in various fields as well. Since the designed algorithm helps to combine the best properties of different criteria for selecting the optimal model, statistical, information-based, residual-based, and others, it overcomes the shortcomings of other approaches, which usually unilaterally consider only one category of criteria. This way, the possibility of genuinely choosing the optimal forecasting model is increased. Consequently, the resulting predictive power leads to noticeably accurate practical results and, hence, the right decisions without too much risk of financial losses. Accordingly, it is believed that all these promising practical implications might justify the slightly higher complexity of the developed mechanism.

CONCLUSION AND FUTURE RESEARCH
In this paper, the new two-stage approach for the ARIMA model selection has been proposed and tested for the real case study of forecasting the road freight transport. In the first stage, the entire family of 16 feasible ARIMA models is identified. In the second stage, the novel FM-MCS mechanism has been applied, which is based on the combination of Mamdani Fuzzy modelling and Monte Carlo simulations. The Fuzzy model was used to access the "hidden" relationships between the ARIMA modelsꞌ measures AIC, HQIC, JB, MASE, MAE, and MSE (Fuzzy inputs) on one side, and the measure %FIT Y = (Fuzzy output) on the other side. The 1000 Monte Carlo simulations were conducted to randomly generate different settings of the fuzzy weights, which has induced the whole cluster of different fuzzy models with different FMs errors for each ARIMA model candidate. The maximal number of 624 occurrences of the minimal FMs errors had occurred for the index * 12 i  , which corresponds to the best model     * * * ARIMA ,1, ARIMA 4,1, 2 12 p q i  . The detected best model provides a promisingly good fit for the real data, particularly considering the major movements of the TS trend. In future work, it is planned to go deeper into research regarding the development of the Fuzzy-based mechanism for model selection, where the other types of the Fuzzy models (e.g., the Sugeno or ANFIS based) are also planned to be engaged. Moreover, the other methods for the optimal setting of fuzzy weights and membership functions are also intended to be deployed (e.g., by applying the Genetic Algorithms). The designed FM-MCS mechanism is also planned to be tested for more complex Box Jenkins and other time series models, e.g., the ones that include exogenous variables (i.e., ARMAX, ARIMAX, SARIMAX models, etc.). Finally, it is planned to test the designed model selection mechanism for the time series that occur in the other areas and/or industries.