Monthly drought prediction based on ensemble models

Non-seasonal model (ARIMA).

Geurts (1977) provide a new forecasting tool known as ARIMA model. It is a modification of the ARMA model, which deals with non-stationary series. According to Mishra & Desai (2005), the general ARIMA model is as follows, (7)${\displaystyle \phi \left(\beta \right){\nabla }^d{y}_t=\theta \left(\beta \right)c_t}$

where, φ(β) is polynomial of Autoregressive (AR) model of order p, θ(β) is polynomial of the moving average model of order q and ∇^d is d^th differencing operator.

Seasonal model (SARIMA).

The SARIMA model is a generalization of the ARIMA model. The major advantage of using the SARIMA model is that it deals with the non-stationary series as well as the seasonal series. According to Durdu (2010) and Mishra & Desai (2005), the general SARIMA model is represented as, (8)${\displaystyle {\phi }_p\left(\beta \right){\varphi }_p\left({\beta }^s\right){\nabla }^d{\nabla }_s^D{y}_t={\theta }_q\left(\beta \right){\Theta }_Q\left({\beta }^s\right)c_t}$

where p and q are the order of non-seasonal part of AR and MA and drepresents the non-seasonal differencing parameter. Besides, P and Q are the order of seasonal part of the AR and MA model, D is a seasonal differencing parameter and srepresents the length of the season.

The development of time series model includes three stages e.g., identification, estimation of parameters and diagnostic check. Initially, the normality assumption of the time series is checked. In case of non-normal series, Box–Cox transformation and log transformation are applied to satisfy the normality condition. After this, Autocorrelation Function (ACF) and Partial Autocorrelation Function (PACF) were used to scrutinize the structure of transformed series and this information is helpful for fitting the appropriate model. So, the initial parameters have been estimated visually from the plot of ACF and PACF of the time series. Then, the parameters of the model have been done by using Maximum Likelihood Estimation (MLE). Various combinations of parameters were used for the fitting of an appropriate model. Among these appropriate models, the choice of a suitable model was based on the minimum value of AIC. Further, the normality condition of residuals was checked by histogram and normal probability plot. We use the R package ‘forecast’ for the time series modelling.

Weighted ensemble drought prediction model algorithm.

The procedure for generating ensemble members in WEDP model is similar to that of the EEDP model Eq. (9). According to Zhang, Dong & Chen (2019), WEDP model algorithm is divided into the following steps:

• The appropriate climate index should be chosen based on the strong correlation with the target month of precipitation (P_i,j,k).

• Assign weights to climate indices through linear a regression model to simplify them into one integrated climate predictor (ICP) φ.

• Calculate the weight of each ensemble member based on the similarity of climate indices between ICP of observed and predicted year.

Traditional weighting procedure.

The impact of climate indices can be signified by an unequal weighting of the ensemble members. Traditional weighting procedure utilizes the Gaussian kernel K(y) to assign a weight to each ensemble member which is based on the similarity of climate indices between the ICP of predicted and observed year (Belayneh, Adamowski & Khalil, 2016; Najafi, Moradkhani & Piechota, 2012; Zhang, Dong & Chen, 2019). A Gaussian kernel function is used to calculate the climate distance between the ICP of predicted and observed year. Therefore, the weight w_i for ith ensemble member is as follows: (10)${\displaystyle w_i=\frac{K\left(\left|{\phi }_i-{\phi }_{n+1}\right|\right)}{\sum _{i=1}^{n}K\left(\left|{\phi }_i-{\phi }_{n+1}\right|\right)}}$

where, φ_i is the ICP of observed year, whereas φ_n+1 is ICP of predicted year, which is obtained from a linear regression model. The Gaussian kernel function is defined as: (11)${\displaystyle K\left(y\right)=\frac{1}{\sqrt{2\pi h}}{e}^{-\frac{y^2}{2h^2}}.}$

The ‘ h’ is the bandwidth kernel parameter and h = σ_ε (σ_ε represents the variance of error between the ICP of observed year and ICP of the predicted year).

Weighted bootstrap resampling procedure.

Smith & Gelfand (1992) proposed the simplest resampling approach known as Weighted Bootstrap Resampling (WBR) procedure. The procedure of WBR is based on random variates via acceptance and rejection method. The procedure of WBR is as follows:

Suppose φ_i i = 1, 2, 3, .., n are samples from g and f. (12)${\displaystyle q_i=\frac{f\left({\phi }_i\right)}{g\left({\phi }_i\right)}}$

and then, (13)${\displaystyle {\omega }_i=\frac{q_i}{\sum _{i=1}^n{q}_i}}$

where, q_i is the ith ratio between f and g. Therefore, w_i, i = 1, 2, …, n are the n weights which are corresponding to ‘ n’ ensemble members and by definition ${\sum }_{i=1}^n{\omega }_i=1$.

Algorithm of conditional ensemble drought prediction model.

According to Zhang, Dong & Chen (2019), the first and second steps of the CEDP model are similar to WEDP model in section 2.3.3.1. The three to five steps of CEDP model are defined as:

• Use the copula function to fit the bivariate distribution between the ICP (φ) and target month of precipitation (P_i,j,k), i.e., C(F_φ(φ), F_Pi,j,k(P_i,j,k)).

• Generate conditional precipitation events 1,000 times given the ICP (φ) of the observed year through the bivariate copula distribution.

• Use the conditional precipitation events i.e., $P_{i,j,k}^1,P_{i,j,k}^2,\dots ,P_{i,j,k}^{1000}$ to generate the ensemble members of the target month.

Copula ensemble.

The basic aim of using the appropriate copula family is to define the dependency structure between ICP and P_i,j,k. The bivariate copula distribution between ICP (φ) and target month of precipitation (P_i,j,k) is as follows: (14)${\displaystyle C\left(F_{\phi }\left(\phi \right),F_{P_{i,j,k}}\left(P_{i,j,k}\right)\right)}$

These, F_φ(φ) and F_Pi,j,k(P_i,j,k) are the marginal distribution of copula function. Marginal distribution has ranged between 0 and 1 because, F_φ(φ) and F_Pi,j,k(P_i,j,k) are the CDF of ICP and P_i,j,k respectively.

According to Zhang, Dong & Chen (2019), the conditional density function is as follows: (15)${\displaystyle f\left(P_{i,j,k}\left|\phi \right.\right)=\frac{f\left(\phi ,P_{i,j,k}\right)}{f\left(\phi \right)}}$

where, f(φ, P_i,j,k) is the joint density distribution between ICP and P_i,j,k. The density f(φ) is the marginal density distribution of ICP. Therefore, the 1,000 times conditional precipitation events are obtained from conditional density Eq. (15), i.e., (16)${\displaystyle P_{i,j,k}^1,P_{i,j,k}^2,\dots ,P_{i,j,k}^{1000}.}$

Finally, the conditional ensemble members are generated by the SPI approach described in section 2.2 on the conditional precipitation events Eq. (16), i.e., (17)${\displaystyle SP{I}_{n+1}=\left(SP{I}_{n+1}^1,SP{I}_{n+1}^2,\dots ,SP{I}_{n+1}^{1000}\right)}$

Copula

In 1959, the general notation of copula was first introduced (Ward, Glorig & Sklar, 1959). A major feature of the Copula function is to examine the dependence structure between the random variables. Following Kao & Govindaraju (2007) and Zhang & Singh (2006), the one-parameter Archimedean copula is expressed as, (18)${\displaystyle C_{\varphi }\left(p,q\right)={\psi }^{-1}\left\{\psi \left(p\right)+\psi \left(q\right)\right\}0<p,q<1}$

The notation ψ⁻¹(⋅) represents the copula generator and it is a convex decreasing function. The notation ϕ is the parameter of Copula function and ϕ(1) = 0. In the present study, four families of Archimedean copulas such as Clayton, Joe, Frank and Gumbel are used.

Accuracy measure tools.

Various tools are used to evaluate the accuracy of predictions (Chen et al., 2014; Pai & Lin, 2005). There are several accuracy measurements tools such as Mean Absolute Error (MAE), Mean Square Error (MSE) and Root Mean Square Error (RMSE), which are scale-dependent measurement tools. Mean Absolute Percentage Error (MAPE) is widely used for measuring the accuracy of predictions of time series models, it is independent of scale. MAPE measures the magnitude of error in term of percentage. According to Parmar & Bhardwaj (2014) and Swanson, Tayman & Bryan (2011), the accuracy measurement tools are described as follows. i.e., (19)${\displaystyle MAE=\frac{\sum _{i=1}^n\left|y_{obs,i}-y_{pred,i}\right|}{n}}$ (20)${\displaystyle MSE=\frac{\sum _{i=1}^n{\left(y_{obs,i}-y_{pred,i}\right)}^2}{n}}$ (21)${\displaystyle RMSE=\sqrt{\frac{\sum _{i=1}^n{\left(y_{obs,i}-y_{pred,i}\right)}^2}{n}}}$ (22)${\displaystyle MAPE=\frac{1}{n}\sum _{1=1}^n\left|\frac{y_{obs,i}-y_{pred,i}}{y_{obs,i}}\right|}$

where, y_obs,i is the observed value, y_pred,i is the predicted value and n is the length of time series.

Normalize Root Mean Square Error (NRMSE) is a non-dimensional form of RMSE and it is the standardized disaggrement between observed and predicted value. It is used to test the predictive accuracy between observed and predicted value. According to Zhang, Dong & Chen (2019), NRMSE is defined: (23)${\displaystyle NRMSE=\frac{RMSE}{y_{obs,max}-y_{obs,min}}}$

where, y_obs,max and y_obs,min are the maximum and minimum value in the observed time series. The smallest value of MAE, MSE, RMSE, MAPE and NRMSE is an indication of better prediction.

Uncertainty measure tool (average bandwidth).

In probabilistic prediction, Chen et al. (2014) presented Average Bandwidth (AB) to measure the prediction uncertainty. The AB provides very important information about the prediction uncertainty., AB is defined as follows: (24)${\displaystyle AB=\frac{\sum _{i=1}^n\left(Q_i^{{}^{1-a∕2}}-Q_i^{a∕2}\right)}{n}}$

where, $Q_i^{{}^{1-a∕2}}$ and $Q_i^{a∕2}$ is an indication of quantiles result of ith the year and a is the level of significance. Therefore, the smallest value of AB is an indication of better prediction.