New Statistical Test to Determine the Best Model for Forecasting Water Discharge Data: Tigris and Euphrates Rivers in Iraq as a Study Case

2.1 Study area

The study area (area of data collection) is located at the water discharge measuring stations located on the Tigris River in the cities of Mosul, Baghdad, and Kut, as well as located on the Euphrates River in the city of Husayabah in western Iraq, as shown on the map in Table 1 and Figure 1.

Table 1. Location of discharge measuring stations

1.png

Figure 1. Area of case study (Discharge measuring stations)

2.2 Time series modeling

Homogeny time series transform to normal distribution (Y_t) is required for stability the variance and improving the normality assumption of the noise series [12].

$Y v, \tau=g \tau(X v, \tau)$              (1)

where: g_τ is the transformation function used to transform it to normal, and then the inverse transformation of Y_t will used to get the generated non-normal series X_t.

$X v, \tau=g-1(Y v, \tau)$                (2)

where: g^-1(Y_v,τ) is the invers transformation function. There are many methods to transformed data to normal, one of these is the original form of Box-Cox transformation [13].

$g(\lambda)=\left\{\begin{array}{cc}\frac{\chi^{\lambda}-1}{\lambda} \\ \ln \chi\end{array} \quad \begin{array}{r}\text { if } \lambda \neq 0 \\ \text { if } \lambda=0\end{array}\right.$                 (3)

in which, value of λ assume from ( -1 to +1, step 0.1) to find optimum λ, it should find the relationship between skewness (cs) and λ with assumption relationship is polynomial second order as shown below:

$\lambda=a_{0}+a_{1} c_{s}+a_{2} c s^{2}$                   (4)

So, the optimum λ where cs = 0, then λ = a₀. To find (a₀, a₁, a₂), it was using polynomial regression between λ_i and c_si by solving these matrices:

$\left[\begin{array}{ccc}n_{\lambda} & \sum C_{S} & \sum C_{S}{ }^{2} \\ \sum C_{S} & \sum C_{S}^{2} & \sum C_{S}{ }^{3} \\ \sum C_{S}{ }^{2} & \sum C_{S}{ }^{3} & \sum C_{S}^{4}\end{array}\right]\left[\begin{array}{l}a_{0} \\ a_{1} \\ a_{2}\end{array}\right]=\left[\begin{array}{c}\sum \lambda_{i} \\ \sum C_{S i} \lambda_{i} \\ \sum C_{S i} \lambda_{i} 2\end{array}\right]$                (5)

2.3 AR models

Autoregressive (AR) models, one of the most important methods in which used for forecasting especially in hydrology and water resources so as to discharge forecasting. It will review the periodic AR models that are modeled with constant coefficients or modeled with periodic coefficients. The periodic normal series (Y v,τ) can be written as [12]:

$Y v, \tau=\mu \tau+\sigma \tau Z v, \tau$        (6)          

where: Y v,τ = gτ (Xv,τ), µτ and  στ are the periodic mean and periodic standard deviation, and (Z v,τ) is the time dependent series with mean zero and variance one.

2.3.1 AR models with constant autoregressive coefficient

The time dependent (Zv_,τ) represented by AR models with constant autoregressive coefficient. For AR (1) model [12]:

$Z_{t}=\Phi_{1} Z_{t-1}+\varepsilon_{t}$            (7)

For AR (2) model:

$Z_{t}=\Phi_{1} Z_{t-1}+\Phi_{2} Z_{t-2}+\varepsilon_{t}$                   (8)

In general, the AR(P) model:

$Z_{t}=\Phi_{1} Z_{t-1}+\cdots \quad \Phi_{P} Z_{t-P}+\varepsilon_{t}$                 (9)

where:

t = (v -1) w + τ

v = no. of years; τ= 1,2, 3…., w; w = number of time interval in years (w = 12 month); Φ_P = the i th partial autocorrelation coefficients; ε_t = the residual or error (assumed normal with mean zero and variance σ²_ε) or (the normal independent variable (white noise) with standard deviation σε). Since (ε_t = σ_ε ². ξ) in which ξ is an independent normal variable with zero mean and variance one [12].

For AR (1):

$\sigma_{\varepsilon}^{2}=\frac{N \cdot \sigma^{2}}{(N-1)}\left\{1-\emptyset_{1}^{2}\right\}$               (10)

$\mathop{\varphi }_{1}={{r}_{1}}\quad $                   (11)    

For AR (2):

$\mathop{\sigma }_{\varepsilon }^{2}=\frac{N.\mathop{\sigma }^{2}.\left( 1+\mathop{\varphi }_{2} \right)}{\left( N-2 \right)\left( 1-\mathop{\varphi }_{2} \right)}\left\{ {{\left( 1-\mathop{\varphi }_{2} \right)}^{2}}-\mathop{\varphi }_{1}^{2} \right\}$                  (12)

$\mathop{\varphi }_{1}=\frac{\mathop{r}_{1}.\left( 1-\mathop{r}_{2} \right)}{\left( 1-\mathop{r}_{1}^{2} \right)}\quad ;\quad \mathop{\varphi }_{2}=\frac{\mathop{r}_{2}^{{}}-\mathop{r}_{1}^{2}}{\left( 1-\mathop{r}_{1}^{2} \right)}$                  (13)

$\mathop{\sigma }_{\mathop{{}}^{{}}}^{2}=\frac{1}{\left( N-1 \right)}\sum\limits_{t=1}^{N}{{{\left( {{Y}_{t}}-\bar{Y} \right)}^{2}}\mathop{{}}^{{}}}$            (14)

2.3.2 AR models with periodic autoregressive coefficient

The time dependent (Zv,τ) represented by AR models with periodic autoregressive coefficients [12].

For AR (1) model:

$Z_{v, \tau}=\Phi_{1, \tau} Z_{v, \tau-1}+\sigma_{\varepsilon, \tau} \xi_{v, \tau}$                  (15)

For AR (2) model:

$Z_{V, \tau}=\Phi_{1, \tau} Z_{V, \tau-1}+\Phi_{2, \tau} Z_{v, \tau-2}+\sigma_{\varepsilon, \tau} \xi_{v, \tau}$               (16)

In general, the AR(P) model as the form:

$Z_{v, \tau}=\Phi_{1, \tau} Z_{v, \tau-1}+\cdots \Phi_{P, \tau} Z_{v, \tau-P}+\sigma_{\varepsilon, \tau} \xi_{v, \tau}$               (17)

where, the $\left(\Phi_{j}, \tau^{\prime} \mathrm{s}\right)$ are the periodic autoregressive coefficient and $\sigma_{\varepsilon}$ periodic variance of the residuals.

For AR (1):

Φ_{1, τ}= ρ _{1, τ;} τ = 1, 2…,12; ρ_{1, τ =} periodic correlation coefficients.

For AR (2):

$\mathop{\varphi }_{1,\tau }=\frac{\mathop{\rho }_{1,\tau }-\mathop{\rho }_{1,\tau -1}\quad.\quad\mathop{\rho }_{2,\tau }}{1-\mathop{\rho }_{1}^{2}}$            (18)

$\mathop{\varphi }_{2,\tau }=\frac{\mathop{\rho }_{2,\tau }-\mathop{\rho }_{1,\tau }.\mathop{\rho }_{1,\tau -1}}{1-\mathop{\rho }_{1,\tau -1}^{2}}$             (19)  

where, τ = 1, 2…,12

The periodic variance of the residual for AR (1) and AR (2) can be written as in equations (20) and (21) respectively [12]:

$\mathop{\sigma }_{\varepsilon \tau }^{2}=1-\mathop{\varphi }_{1,\tau }.\mathop{\rho }_{1,\tau }$                  (20)

$\mathop{\sigma }_{\varepsilon \tau }^{2}=1-\mathop{\varphi }_{1,\tau }.\mathop{\rho }_{1,\tau }-\mathop{\varphi }_{2}.\mathop{\rho }_{2,\tau }$           (21)

4.1 Results of data generation

To generate the monthly flow discharge, a computer programme was written, using the available data for the gauges (Mosul, Baghdad, Kut, and Husaybah). The purpose is to generate a sample of 50 years using the AR (1) and AR (2) models and with two cases of a constant autoregressive coefficient and a periodic autoregressive coefficient for each model, and for the four gages. The programme is written in the form of a subroutine. The first subroutine includes a test of the homogeneity of the data. It includes checking the trend existing by using t-test Eq. (23) and t-table with a significant level of 0.975 and degree of freedom equal to (n-2) for the four gauges is 1.96, and comparing it with the t-calculated. It was found no trend exists for the four gauges because the t-calculated is less than the t-table, so there is no need to remove the trend and the data is homogeneous for all gages.

The second subroutine consisted of the test of normality of data. It was found the skewness (cs) for all gauges was greater than zero as shown in Table 2. Therefore, there was a need to transform the data using Eq. (3) to evaluate the optimum (λ) after finding the coefficients of the second-order polynomial equation for all gauges through Eq. (4) and Eq. (5). The optimum (λ) for the gauges is illustrated in the Table 2.

Table 2. Skewness and optimum (λ) for all gauges

The third subroutine of the program, the periodicity, was analysed for six harmonics to determine µτ and στ using Eq. (24) to Eq. (27), using the variance analysis and F-test to indicate whether the periodicity exists or not. If the periodicity exists, we will find it µτ and στ, after specifying the harmonic number that is significant through explaining percent, we will remove the periodicity Eq. (28).

The test mean and standard deviation of the results should be zero and one, respectively, or they will be standardized. The result is (ε), which is the residuals after removing the periodicity and standardizing. From variance analysis, it was found that the periodicity existed for all gauges except Baghdad gauge, which didn’t have periodicity, and according to explain percent, it was found that the number of harmonics significant was 1,2,5, and 4. For example, the harmonics significant were 1,2,4, and 3. In Table 3 are also shown Fourier coefficients of fitted harmonics to means and standard deviations.

Finding the autocorrelation for the residual (ε). This was done in the fourth subroutine of the program, from lag equal one to lag (n/4). Figure 3 illustrates the relationship between lags and r_k (ε).

3.png

Figure 3. Relationship between r_k (ε) l and lags

From Eq. (11) and Eq. (13) it was determined the partial autocorrelation function for AR (1) and AR (2) models as constant. The results are shown in Table 4 for all gauges.

The fifth subroutine of the program, including the calculation of the periodic autocorrelation function with lag 1 and lag 2 for all gages, is illustrated in Table 5, which leads to the computation of the partial autocorrelation function for AR (1) and AR (2) models as periodic Ф. The results are shown in Table 6.

Table 3. Test periodicity, no. of harmonics, and Fourier coefficients for means and standard deviations

Table 4. Constant Ф1 for AR (1), constant Ф1 and Ф2 for AR (2)

In Tables 7 and Table 8 it was calculated the variance (σε ²⁾ of dependent stochastic component using AR (1) and AR (2) as constant Ф by Eq. (10) and Eq. (12), regarding periodic Ф the calculation of (σε ²) was used Eq. (20) and Eq. (21), this is done in sixth subroutine of the program.

In subroutine seventh through a tenth of the program and from Eq. (9). It was computed (ε) for AR (1) and AR (2) models with constant Ф, and periodic Ф has used Eq. (17).

Using 600 standard normal random number for generating 50 years with eliminating 5 years data to avoid an error that produce due to an assumption of initial (z). It is done by Eq. (9) and Eq. (17) with substituting (εt = σε2. ξ) and (εt τ = σε τ 2. ξv,τ) respectively. Figure 4 to Figure 7 represented the data generated for all gages and for AR (1) and AR (2) models with constant Ф and periodic Ф. Discharges were generated for Mosul gage up to the year 2042, for Baghdad gage up to the year 2036, for Kut gage up to the year 2047, and Husaybah gauge up to the year 2038.

4.png

Figure 4. Data generation for Mosul gauge

5.png

Figure 5. Data generation for Baghdad gauge

Table 5. Periodic autocorrelation function

Table 6. Partial autocorrelation function

Table 7. Values of σ_ε² for AR models with constant Ф for all gauges

6.png

Figure 6. Data generation for Kut gauge

4.2 Results of a suggested new test

Previously, it was looking for a method to choose the best model. Regarding the portmanteau lack of fit test, the result from Eq. (31) is illustrated in Table 9 regarding the results for Mosul gauge. Although the Q-calculated is greater than the chi-square Table, that means all models are not accepted. Therefore, it should be attempting to test high-order models such as AR (3). However, it could choose AR (2) with periodic Ф as a better model because it has a lower Q-Calculated. Also, from the Table 9, it is found that the AR (2) with periodic Ф is suitable for Baghdad. AR (1) with periodic Ф is a better model for Kut gauge according to lower Q-calculated.  Also, Husaybah gauges could be decided AR (2) with periodic Ф. The best model was chosen because the Q-calculated and chi-square Tables were the lowest, although all models were suitable. Regarding the AIC test, the results of Eq. (30) are shown in Table 10 with the note that the mean of (σε2) when calculating the AIC for AR (1) and AR (2) with periodic, Ф, the AR (2) with periodic Ф model was better for Mosul, Baghdad, and Kut gauges. According to the lower AIC values, Husaybah's AR (1, with constant) was the best.

7.png

Figure 7. Data generation for Husaybah gauge

Table 8. Values of σε 2 for AR models with periodic Ф for all gages

Table 9. Port Manteau Lack of fit test

Table 10. AIC test results

8.png

Figure 8. Autocorrelation r_k(ξ) of AR (1) with constant Ф, for Baghdad gauge

9.png

Figure 9. Autocorrelation r_k(ξ) of AR (1) with periodic Ф, for Baghdad gauge

10.png

Figure 10. Autocorrelation r_k(ξ) of AR (2) with constant Ф, for Baghdad gauge

Suggestion test depending on the rk (ξ) values, the drawing between rk (ξ) and lags (k) produces the correlograms. Correlograms for all models and for all gauges are shown in Figures 8 to 14 with confidence limits that result from Eq. (29). The objective is to find the number of points that lie outside the confidence limits and could be decided upon the better model according to the lower points that lie outside the limits. Firstly, it will compute the points outside the limits using Eqns. (32) to (35) and compare with Correlogram graphics will also be compared with the results from the Portmanteau Lack Test and the results from the AIC Test for all models and for all gages. Through Correlograms from Figures 8 to 11 for the Baghdad gauge, it is very clear that the number of points that lie outside the limits was 1,2,1, and zero for the autocorrelation function of residual of AR1 with constant coefficients, AR1 with constant, and AR2 with periodic, Ф respectively.

11.png

Figure 11. Autocorrelation of AR (2) with periodic Ф, for Baghdad gauge

12.png

Figure 12. Autocorrelation r_k(ξ) for Mousl gauge

From Eq. (32) through Eq. (35) will compute the number of points that are outside of the limits. The result is illustrated in Table 11, where the AR (2) with periodic Ф is the better model and that is identified with the graphics and also the same conclusion when using the portmanteau lack test and when using the AIC test. Correlograms for Mosul, Kut, and Husaybah were shown in Figures 12 through 14. The results for all gauges when using the suggested test are illustrated in Table 12.

It was seen that there were no suitable models for Mosul gauges because the AC% was less than 95% for all models, so we should try analysis higher-order, anyway. Regarding the AR (2) with periodic Ф model is accepted, this identical with portmanteau lack test and when used AIC test. The better model for Kut gauge is for AR2 with constant. This is identical to the portmanteau lack test and doesn’t match when using the AIC test. So, the suggested test is more reliable if compared with the graphics. Regarding Husaybah gauge from the suggested test, the better model is AR (2) with periodicity. This is Cleary from graphics and similar when used in portmanteau tests and far away when used in AIC tests. Hence, the AIC test gives the AR (1) with constant So the results that are produced from the suggested test are reliable when compared with graphics. This is the proposed statistical method for choosing the best model that depends on the value of the autocorrelation coefficient of the residual rk(ξ) of lag through the acceptance percentage, which represents the percentage of the number of points within the upper and lower limit.

13.png

Figure 13. Autocorrelation r_k(ξ) for Kut gauge

Table 11. AIC test results

Table 12. Results of suggested test

Table 13. Results of suggested test

14.png

Figure 14. Autocorrelation r_k(ξ) for Husaybah gauge