A hybrid forecasting model based on the group method of data handling and wavelet decomposition for monthly rivers streamflow data sets

2.1 The WD

The WD algorithm combines the traditional AI forecasting algorithms developed in past periods. A novel methodology to innovate the TS forecast model and the WD algorithm has emerged as an ever more useful and vital tool for the models [33]. The WD established a predicting method for TS, which presented to develop several decomposition goals that utilize signals for repetitive trous wavelet transformation and benefit from being a shift-invariant. Further, the WD algorithm can de-noise nonstationary signals into sub-signals at different levels with a relevant resource for improved streamflow elucidation [13,27]. The mathematical structure of the WD is described in Eq. (1) as a continuous TS x(t) (∵t ∈ [−∞, ∞]).

where t is the time, τ  is the time step that iterates the window function,  and s ∈ [0, ∞] for wavelet scale defined, respectively. As ψ(t) is known as the mother wavelet and expressed in Eq. (2),

whereas the continuous wavelet transform (CWT) illustrated in Eq. (3) [28]

where ψ ¯ ( t ) and W(τ, s) described to complex conjugation and consumed the total time of the inputs multiplied by the scale and shifted the form of the wavelet function ψ(t). Forecast by the CWT is non-feasible since wavelet-coefficient estimation at each credible scale is time-consuming and provides some outcomes. Thus, the WD is preferred in most prognostic models due to the time-efficient and simple algorithm. The WD consists of selecting the scales and position concerning the powers of 2, alleged dyadic scales, translations, and the interpretation that may be most precise and effectual. The significant impact of implying the WD is its robust algorithm because it has no hypothetically specious or parametric testing technique [34]. The WD is described as a mathematical structure in Eq. (4)

where m, n ∈ Z switching as the scale and length of the TS, respectively, s ₀ > 1 is defined for fixed dilation step, and τ ₀ > 1 is the position parameter. The parameters s ₀ = 2  and  τ ₀ = 1 are the most common preference. For a discrete TS x(t) that occurs at t discrete time, then WD formation is defined as follows:

In Eq. (5), x(t) is defined as discrete TS(∵1 ≤ t ≤ N − 1) on the power of 2 (N = 2 ^m ), and the wavelet coefficient of the WD (W _m,n ) is described on s = 2 ^m   and  τ = 2 ^m n scales, where n is the time-translation component that varies within the ranges 0 < n < 2M − m  and  1 < m < M, and M is the estimation of decomposition levels, which is interpreted the TS into numerous wavelet input segments on the particular decomposition level.

Conclusively, the desired decomposition level of the time series in WD plays a significant factor in retaining figures and decreasing the time series curvature. Though, there has no such hypothesis regarding the number of decomposition levels for data sets. Formation (6) is used to determine the number of decomposition levels [34].

According to Mallat's theory, the discrete TS x(t) decomposed into a new series of linearity independent estimate and detail signals through the inverse WD, which is described in Eq. (7) [13,35];

or defined as

Eq. (8) A _M (t) is the residual term at M decomposition level and D _m (t) has detailed subseries (∵m = 1,2, …, M) that are taken as small features in interpreted value data. The following various forms of TFs have emerged and are vastly used in the model [36]; the efficiency of the conventional GMDH models can be improved.

Three different TF types (9), the sigmoid, radial basis, and hyperbolic tangent functions have been proposed and improved the conventional GMDH model's performance and proved the numerous substitutes into a single network. In conventional GMDH models, most researchers have measured the partial quadratic polynomials (PQP) as the TF and demonstrated that utilizing various TF exclusively in a model gives preferred outcomes rather than using the units for a single TF [37,38,39]. Nevertheless, multiple TF (9) are accessible for the modified GMDH model, and it has provided more reliable and accurate outcomes than the single TF models.

2.2 GMDH

The GMDH model has been established to display the complex nonlinear framework by utilizing a TF to communicate the connection between input and output data sets. The Volterra linear function (10) describes the input and output mapping with a neural network like the GMDH model.

The mathematical structure (11) comprising the two order variables (neurons) is expressed in the following form:

In mathematical structure (11) have six coefficients {b _o , b ₁, b ₂, b ₃, b ₄, b ₅} in the g function, and these are unknowns [37]. Let us consider the coefficients are determined to expend the least square method. The system's observed input and forecast (output) variables are set for y and z, respectively.

The following are the iterative steps of the MGMDH model structure

Step I: Initially considered M ₁ = m as the neurons of the input vector in the first layer Y = (y _i , y _i , …, y _n ), whereas m is the number of inputs and taken j = 1 as a threshold.

Step II: Set N = M j ( M j − 1 ) 2 for all independent variables Y as new variables x ₁, x ₂, …, x _N in the data set and by the hypothesis of the PQPs constructed Eq. (10).

Step III: The aggregate number of partial descriptions (PDs) at the current layer has equivalent to the following combination:

Eq. (12) n is the number of nodes that have been chosen input variables from the nodes of the previous layer. The sort of PDs in each layer of the network is similar.

Step IV: Estimated the coefficients of the TF by the least squared (MSE) in the form of

where B = {b ₀, b ₁, …,  b ₅} has the input vector of the coefficients of the variables and Y = {y ₁, y ₂, …,  y _M } ^T is the output value from the following observation (14)

Step V: Distinguish the absolute best neuron out of T neurons. If the tiniest estimation by RMSE is less than the threshold, then the input will be terminated, or new input variables are introduced {y ₁, y ₂, …, y _n , y′} (since y′ specified for the RMSE estimate) and T = T + 1 will be set; otherwise, the GMDH model Steps II and III will be repeated.

Step VI: The termination standard for the iteration of the model is as follows:

The lowest estimation at the individual layer acquired this iteration is associated with the lowest values found at the previous level; the lowest value is compared to the smallest assessments using the GMDH-constructed layer model. Unless the improvement is accomplished, till that repeats and rehashes steps II–IV. The last layer estimated the one node classified by the best execution, selected as the output node, and rejected the rest of the layer nodes. Finally, the GMDH model forecasting estimation is acquired, and then, the iteration is terminated and acknowledged by the network.

2.3 The WGMDH model

The methodology for a hybrid forecast model using the GMDH model combined with WD is anticipating future data sets by utilizing wavelet capability as a time-frequency restriction. Consequently, the best possible data sets of the GMDH algorithm have been developed through decomposition with a few good contrast strategies in this forecasting algorithm. The trous WD lets a significant translation of the filter outputs. The A-trous wavelet has been shown and begins with performing successive convolutions with the discrete low-pass channel.

whereas D ₀(t) = Y(I) has the original time series. The development of separations between instances has priority and 2 ⁱ I illustrated that a trous name is synonymous with this algorithm, and the low-pass filter is a spline. This expression has a point symmetric; it has also helped to minimize; however, lastly, it does not consider how time is essentially asymmetric. The smoothed structure and contrast of progressive variants enable the w _i wavelet coefficients to be accomplished by the detail signal.

Eq. (16) is contained detailed signal or wavelet coefficients with a small detention feature of interpretational in the data value. The wavelet coefficients are iterated into the GMDH and analyzed as inputs. The WGMDH model ability was carried out by choosing the right inputs to overcome by the hybrid algorithm for the appropriate estimations.

2.4 Statistical performance

The use of statistical estimations is crucial for the reliability and credibility of any forecasting algorithm. Therefore, the estimations of the forecasting algorithms have been investigated through three statistical methods, such as the Mean Absolute Error (MAE), the Root Mean Square Error (RMSE), and the correlation coefficients (R). The statistical estimation outcomes fitted to the best fit line and the forecast performances estimated by the following indices [38]

where y t , y ˆ t , y ¯ t and y ¯ ˆ t are the actual (streamflow observed), forecasted, average streamflow observed, and average forecast values at a time t and, respectively, whereas m is the total number of observations.

The normalization and performance of the data set are used in the training and testing process and observed the estimations to select the best model with the least values of MAE and RMSE and the highest value of R (effect of flow).

3.1 Indus River

The Indus River, locally called Sindhu discovered by the Scylax of Caryanda, the Greek advisor of the King Darius of Persia in 515 BC, and the early classical age information of the river had spread out in the West. The Indus River is one of the largest rivers in Pakistan and the Asian subcontinent and plays a significant role in the agricultural production industry. This river flows by Western Tibet, China, and India along with the full length of Pakistan to the south direction (via Gilgit-Baltistan) and falls into the Arabian Sea near the port city of Karachi in Sindh.

Figure 2

Indus river location of the study area.

The Indus River streamflow 484 months rate has been collected from the Department of Irrigation and Drainage, Sindh, from September 1978 to December 2018 (Figures 2 and 3).

Figure 3

Streamflow of Indus River (484 months rate from 09-1978 to 12-2018).

3.2 Chenab River

As five main rivers are flowing in the Punjab province of Pakistan, one of them is the Chenab River which originates in the upper Himalayas of the Indian state of Himachal Pradesh. It is flowing from Jammu and Kashmir territory into the country's fields before flowing into the Indus River near the city of Uch Sharif. Water has been granted for this significant River under the Indus Water Treaty.

Figure 4

Chenab river location of the study area.

The Chenab River streamflow 550 months rate has been collected from the Department of Irrigation and Drainage, Punjab Barrage from March 1973 to December 2018.

Figure 5

Streamflow of Chenab River (550 months rate from 03-1973 to 12-2018).

The hydrological data sets are used to forecast the river streamflow of both well-known rivers and validate the modified GMDH algorithm for mean monthly river streamflow; then, the same data sets are used to develop hybridized wavelet decomposition GMDH model (Figures 4 and 5).

3.3 Model structure

Appropriate picking up the streamflow TS input has one of the most significant steps in constructing a new forecast prototype. In case of the average discrepancy ratio either less (underestimated) or more than (overestimated) the unit, the following linear transformation (LT) (20) estimated the inputs for the model with decomposed components which are scaled at (0, 1)

In LT (19) has been represented y _t is the initial streamflow data set decomposed by Empirical Mode Decomposition (EMD) with h _t = X _t − min(X) the normalization of the data is conducted before all the Intrinsic Mode Functions (IMFs) are observed by the Encoder-Decoder Long Short Term Memory (En-De-LSTM) model at the time t and X _t the decomposed components and h _max = max(X) − min(X) the maximum values between the actual values.

The corresponding input variables allow the network to get the aspired output and prevent essential information trial effectively. There are no fixed-designed rules for the TS input variables to develop the forecast prototype, while a general structure can be followed based on the successful implementation of water resources issues [39,40,41]. In this research work, six input model combinations have been developed to explore the input-output variables model, equivalent to lagging variables from the previous monthly streamflow time series cycles y _t−p , where 1–6 months are set for p, respectively. Eqs. (21) and (22) mathematically model structures followed for original y _t−p and DW _t−p monthly streamflow statistics.

and

where y _t−1 corresponds to the last month's streamflow value, there DW _t−1 effectively estimated the components D _s (D ₂, D ₄, and D ₈) by evaluating the y _t−1 and the input model structure values described in Table 1 for monthly streamflow forecasting of the Indus River and Chenab River (Figures 6a and 6b).

Figure 6

(a) Flow structure of the study; (b) flowchart of proposed WMGMDH model.

4.1 Fitting GMDH model to the data sets

It is significant to pick the proper input series for the GMDH model and then use it; therefore, the six(M1−M6) input structures are suitable for the model used for streamflow data sets in two-tier training and testing phases. Significant and various studies have shown the positive results of the radial basis function (RBF) in the use of the AI models for the configuration and estimation of TS data sets [16,29,30]. Therefore, this study uses the RBF kernel function to control the streamflow modeling and TS forecast system. There are two parameters (γ, σ ²) used for the estimation of the GMDH model. These parameters are used to test the model more efficiently and systematically to address possible shortcomings in the estimation and trails method by the standard grid search procedure. The model performance has been analyzed to determine the parameter values by grid exploring γ  and σ ² in the range of 10–1,000 and 0.01–1.0, respectively, and the cross-validation procedure is applied for calibrating parameters to evade over-fit. The k-fold cross-validation of the training phase for approximation forecasting error is given for each pair of hyper-parameter (γ, σ ²) in the search space. The forecasting best-fit algorithm is incorporated based on execution and performance standards.

The R, MAE, and RMSE statistical estimates for the training and testing phases from a combination of six input specimens (M1−M6) of the GMDH model are presented in Table 2 for the Indus and Chenab rivers. It is observed that the training and testing phases M6 and M5 input models R had the greatest estimate whereas the least estimations have for MAE and MSE; the testing estimations of M6 and M5 input combinations of R (0.9635 and 0.9547), MAE (0.0247 and 0.0290) and RMSE (0.0011 and 0.0015) shows better estimations for the model.

4.2 Fitting hybrid model wavelet GMDH to the data sets (WGMDH)

A new forecast streamflow algorithm has been developed with a crossbred thought that uses the wavelet and GMDH algorithms. The primary TS data sets have been decomposed into WD periodic components using Mallat's DWT formation before applying the GMDH algorithm. As chosen, the decomposition levels based on the perceived time series have been decomposed into wavelet components. The best decomposition of TS data sets in wavelet analysis plays a significant role in conserving and minimizing information distortion.

The number of decomposition levels is estimated at M = 2.6848  and  M = 2.7403 for the Indus and Chenab rivers using n = 484  and  n = 550 monthly streamflow in the formula (6). Still, the data sets decomposed on round-off value M = 3 (i.e. 2–4–8 months) decomposition levels used for both rivers estimations. The efficacy of WDs is measured by the correlation coefficients and analyzed between streamflow TS data sets and the various decomposition levels of wavelet coefficients. The measurement of the correlation coefficients of each TS WD and primary monthly streamflow data are described in Table 3. It is observed that the D ₁₁ (2-months mode) has a low correlation coefficients value compared with  D ₂₁ (4-month mode) and  D ₃₁ (8-month mode) components, whereas these have well shown higher correlations. Subsequently, the correlated WDs have been well shown, such as D _2p ,  D _3p ,  and  A _3p (approximate mode) components added to each other for developing the new input series and that developed series used with the WGMDH model as GMDH inputs. The proposed model structure of WGMDH is presented in Figure 7, and Figure 8 shows the combinations of useful details and approximation components mode such as D _1p ,  D _2p ,  D _3p , A _3p , and (D _2p +  D _3p + A _3p ). The six variant specimens of the newly input series data described in Table 1 are used to forecast as in the previous GMDH model.

Figure 7

Framework of the study.

Figure 8

Wavelet decomposed combinations of effective details and approximation components mode of TS data sets of Indus and Chenab Rivers.

The computational codes for GMDH, WGMDH, LSSVM, and ANN models are written in the MATLAB software with the wavelet toolbox, which communicates the streamflow forecast and models estimations are presented in Tables 4–7 in the training and testing phases. The common critical approach concerning the machine learning model's achievement is the training and testing phases. The training phase is executed to the enhanced characteristic of the developed model algorithm while the testing phase estimations are used for the model validation and once the model is validated, then the training phase is eliminated from the testing phase [42,43,44]. The best statistical estimates are selected from the testing phase to appropriate forecast analysis; therefore, the testing phase estimations have the best forecast outcomes for developing the model. Further, the WGMDH model reveals the strong effects of the Indus and Chenab Rivers streamflow forecast, clearly described in Table 4. In view of the table analysis, the input model WM5 has estimated the most considerable and best for the Indus River; the largest and lowest values are R = 0.988 and MAE = 0.0117, RMSE = 0.0004 in the training phase, while the same input model has been observed that the highest and least optimal estimations as under R = 0.9859 and MAE = 0.0122, RMSE = 0.0002 in the testing phase. Similarly, in the table, the input combination WM4 has given the best estimate for Chenab River, with the largest and lowest values R = 0.9901 and MAE = 0.0119, RMSE = 0.0005 in the training phase, while the same input model has been opted for the highest and least estimations as under in the testing phase R = 0.9798 and MAE = 0.0127, RMSE = 0.0008.

4.3 Comparisons of forecasting models

The developed WGMDH model has been compared with GMDH, LSSVM, and ANN models based on performance and accuracy inspection by statistical estimations in the testing phase. The testing phase of Table 7 is reported that the least statistical estimations of MAE and RMSE have in the WGMDH model compared to other models; this is evidence of the valuable performance. Furthermore, the same model also provided authoritative evidence for estimating the most considerable value R; therefore, the WGMDH model is significant and better than the other models.

The correlation coefficient (R) estimated on Indus and Chenab Rivers for the GMDH model is 0.9635 and 0.9547, for LSSVM model is 0.9335 and 0.9467, for ANN model is 0.9125 and 0.9144 and for WGMDH model is 0.9859 and 0.9798 highest correlation value estimates and shows averagely 2.24% and 2.5% increased, respectively. It is a good sign for the developed model with the combination of WD and GMDH algorithms as compared to the conventional models in terms of model accuracy. In addition, the MAE and RMSE on Indus and Chenab Rivers data sets are statistically the smallest values in the forecasting models and Table 7 illustrated that the developed model had achieved the smallest value compared to the other models; this is again a valuable element for the model efficiency.

Figures 9 and 10 are demonstrated the predicted linear trend line scatter plotting for the WGMDH, GMDH, LSSVM, and ANN models from the testing phase and clearly is described the outperformance as relatively close to observed data on both streamflow stations by the WGMDH model. Therefore, the application of WGMDH in river streamflow forecast is better than that of conventional GMDH and other models. The scatter graphs are described by the regression fit line equation y = a + bx using MATLAB application, where a and b are the coefficients and defined for the GMDH, WGMDH, LSSVM, and ANN models, respectively. The WGMDH model has fewer scattered estimates, and the R-value of the WGMDH model has close 100 per cent (i.e. R = 0.9895 and R = 0.9798) for both Indus and Chenab data sets compared to the GMDH, WLSSVM, WANN and WMARS models. Therefore, this study concludes that the WGMDH model is superior to other models as it provides better forecast estimations for both streamflow stations.

Figure 9

Predicted linear trend line scatter graphs for the Indus streamflow by WGMDH, GMDH, LSSVM and ANN models during the testing phase.

Figure 10

Predicted linear trend line scatter graphs for the Chenab streamflow by WGMDH, GMDH, LSSVM, and ANN models during the testing phase.