This study introduces ensemble empirical mode decomposition (EEMD) coupled with the autoregressive integrated moving average (ARIMA) model for drought prediction. In the realm of drought forecasting, we assess the EEMD-ARIMA model against the traditional ARIMA approach, using monthly precipitation data from January 1970 to December 2019 in Herat province, Afghanistan. Our evaluation spans various timescales of standardized precipitation index (SPI) 3, SPI 6, SPI 9, and SPI 12. Statistical indicators like root-mean-square error, mean absolute error (MAE), mean absolute percentage error (MAPE), and R2 are employed. To comprehend data features thoroughly, each SPI series initially computed from the original monthly precipitation time series. Subsequently, each SPI undergoes decomposition using EEMD, resulting in intrinsic mode functions (IMFs) and one residual series. The next step involves forecasting each IMF component and residual using the corresponding ARIMA model. To create an ensemble forecast for the initial SPI series, the predicted outcomes of the modeled IMFs and residual series are finally added. Results indicate that EEMD-ARIMA significantly enhances drought forecasting accuracy compared to conventional ARIMA model.

  • Improved drought forecasting: Our study introduces the ensemble empirical mode decomposition–autoregressive integrated moving average (ARIMA) model, enhancing drought forecasting accuracy over traditional ARIMA methods.

  • New drought model for Western Afghanistan: We present a customized model for drought prediction in Western Afghanistan, based on the standardized precipitation index.

Drought, which is generally known as a shortage of rainfall below the average for an extended period of time, is one of nature's most catastrophic disasters. It leads to adverse influences on both the environment and human life. It is vital to accurately predict the quantity of rainfall and drought because it can negatively affect agricultural products and water reservoirs due to features of high frequency, prolong occurrences, widespread influence, and detrimental consequences (Khan et al. 2020; Yang et al. 2020). To achieve an accurate prediction of drought in a region, precise rainfall prognostication seems obligatory. Drought research can benefit from more precise and timely rainfall predictions, which can also greatly enhance future water management policy in many ways. Due to its location in a dry region, Afghanistan must pay close attention to this issue. For instance, Afghanistan endured its most destructive droughts in the previous 50 years. Major cities including capital Kabul experienced critical water shortages during past decades, thousands of villages ran out of drinking water, surface water flow fell down significantly, and dams and reservoirs were forced to operate at their bare minimum water transport capacity due to low inflow and high temperatures (Jadin 2018; International Federation of Red Cross and Red Crescent Societies 2021).
Figure 1

Location of precipitation stations, Herat province, Afghanistan (the study area).

Figure 1

Location of precipitation stations, Herat province, Afghanistan (the study area).

Close modal

The gradual onset and development of drought can be beneficial for researchers to predict these occurrences ahead of time (Rossi 2003; Cancelliere et al. 2007). Mitigating risks has become an essential aspect of planning and policy-making processes. In recent years, this has been facilitated by the emergence of various modeling techniques that provide valuable information about precipitation shortfalls and ultimately enhance drought monitoring capabilities. Because drought adversely affects crops seasons and exacerbate water shortage crisis and have direct influence on human quality life, an early authentic drought prediction plays an indispensable role to mitigate its side effects (Alam et al. 2014).

Through investigating of drought literature, it can be easily found that there are several drought indices (DIs) including standardized precipitation evapotranspiration index, Palmer drought severity index, standardized precipitation index (SPI), effective drought index, and reconnaissance drought index. SPI is one of the most common DIs widely employed for drought forecasting (McKee et al. 1993) and has been endorsed by the World Meteorological Organization (WMO) (Hayes et al. 2011). The calculation of SPI is based on distribution of precipitation, with monthly rainfall amounts being normalized using a gamma or Pearson type III distribution (McKee et al. 1993). Several benefits for employing the SPI index including less complex computation process, offering stable spatial interoperation in different climate area, and its probabilistic features (Guttman 1998; Zargar et al. 2011) were among compelling reasons to utilize this index in the current study. The importance of SPI is mostly in areas with data-sparse characteristics where immediate access to other parameters including evaporation, streamflow, and soil moisture information is impossible (Pasteris et al. 2005). In addition to the aforementioned advantages of using SPI, it can be computed in different timescales that offers a wide perspective into various kinds of drought. For instance, short to medium timescale proves effective in understanding meteorological and agricultural drought, while longer timescales are more applicable to hydrological and socioeconomic drought (Gumus & Algin 2017). The SPI possesses a multitude of versatile applications. In recent years, it has been employed to evaluate groundwater drought and to perform spatiotemporal analysis of floods and droughts (Kumar et al. 2016; Liu et al. 2018; Anshuka et al. 2019).

Although forecasting drought events present inherent challenges due to the complexity of time series data, several approaches have emerged in recent years, including autoregressive integrated moving average (ARIMA), artificial neural network (ANN), support-vector machine (SVM), fuzzy logic, and Gaussian process regression.

For instance, a singular spectrum analysis (SSA) and single least-square support-vector machine (LSSVM) model combined to anticipate drought. Data on the SPI are obtained for the southern Taiwanese Tseng-Wen reservoir catchment (1975–2015). According to the results, LSSVM1 (an LSSVM-based technique that uses antecedent SPI as an input) is not as effective as SSA-LSSVM2 (an LSSVM-based model combined with SSA utilizing the antecedent SPI as input). Furthermore, SSA-LSSVM3 is the most appropriate model for SPI drought forecasting when SSA-LSSVM2 and SSA-LSSVM3 (an SSA-LSSVM-based technique that uses antecedent cumulative monthly rainfall as an input) are compared for performance (Pham et al. 2021; Alawsi et al. 2022). In addition, the wavelet packet transform (WPT) method was applied to the data preprocessing step. The model's construction and evaluation were based on India's SPI from 1985 to 2013. The results show that hybrid strategies – WPT-ANN and WPT-support-vector regression (SVR) – perform better than ANN and SVR, which are standalone machine learning approaches. The hybrid WPT-ANN model fared better than the WPT-SVR model for the majority of grid sites (Das et al. 2020). In Ethiopia's Awash River Basin, three models were utilized to forecast SPI 3 and SPI 12 between 1970 and 2005: SVR, ANN, and hybrid wavelet and ANN (W-ANN). The performances of each model were evaluated using R2, root-mean-square error (RMSE), and mean absolute error (MAE). The findings show that the most effective method for predicting droughts is to use wavelet neural networks (W-ANN) (Belayneh & Adamowski 2012).

Among these methods, ARIMA, a data-driven approach, has been widely utilized in drought and streamflow forecasting due to its ability to consider serial correlations inherent in time series data (Mishra & Desai 2005). Compared to alternatives like ANN and exponential smoothing, ARIMA offers advantages such as a comprehensive search stage involving model identification, estimation, and diagnostic checking, making it a popular choice for various types of time series data (Zhang 2003; Mishra & Desai 2005). However, it is important to note that according to Nourani et al. (2009), ARIMA is most suitable for stationary and some simple nonstationary datasets, and it may have limitations in capturing nonstationary time-oriented data. To address this limitation, decomposing techniques have been introduced to enhance the accuracy of ARIMA models for nonstationary series. For example, in a study by Rezaiy & Shabri (2023), a hybrid wavelet and ARIMA (W-ARIMA) model was employed for drought prediction in Kabul, Afghanistan. This approach aimed to overcome the limitations of the ARIMA model in capturing the nonstationarity of SPI values. The study's results demonstrated the superior performance of the W-ARIMA model over the traditional ARIMA in predicting drought events, as evidenced by various error metrics. Similarly, another study (Shabri 2014) focused on utilizing the SPI time series as a drought indicator. The investigation assessed the effectiveness of a coupled wavelet-adaptive neuro-fuzzy inference system (WANFIS) model for drought forecasting. The analysis incorporated SPI 3, 6, 9, and 12 series to evaluate drought conditions. To illustrate the applicability of the WANFIS model in drought forecasting, rainfall data from the Klang River basin in Malaysia were obtained. Time series data typically include both linear and nonlinear fluctuations from many applications. Low-precision linear ARIMA and nonlinear ANN models can be used to model discrete time series of drought events. But when the advantages of ARIMA and ANN models were merged, a hybrid model was created that performed better than either model alone. In the study by Khan et al. (2020), a novel ARIMA-ANN hybrid forecasting model is created, utilizing the discrete wavelet transform. It is determined that for SPI and standard index of annual precipitation, which are regarded as meteorological DIs, the wavelet-ARIMA-ANN (W-2A) model assisted in enhancing forecasting accuracy. In the Algerois basin of Algeria, the SPI at three timescales (SPI 3, SPI 6, and SPI 12) was projected using the combined wavelet (W) and neural networks (W-ANN) technique (1936–2008). The solo, ANN, and two conventional stochastic models (ARIMA and SARIMA) were contrasted with the combined model. According to the findings, the wavelet approach enhances the quality of the data, and W-ANN performs best when considering various statistical parameters (Djerbouai & Souag-Gamane 2016).

These techniques involve generating stationary sub-series that are predicted using ARIMA, and the reconstructed combination of these sub-series serves as the prediction for the original time series. Several studies showed that hybrid models are more capable than applying individual models in time series analyzing except when the single model decomposed by preprocessing techniques such as wavelet transform. However, the application of different hybrid models increased in time series analysis because of their higher prediction accuracy opposed to single models, and some research studies announced that extreme utilization of such hybrid models lowered prediction accuracy and produced a model with no significant excellency (Khan et al. 2020).

Hydrological time series records including drought events encompass some subprocess and modeling, making it occasionally improper to use an individual model. To achieve a better comprehension of hidden physical features in time series, the original data's time–frequency change laws can be investigated by applying decomposition techniques (Wang et al. 2018). To simplify the prediction of a complicated series, a suitable decomposition approach can be exercised to the original time series to separate a single series into some sub-series, capturing a different aspect of original signal, and then to improve precision of the forecasting accuracy.

A majority of decomposition techniques used in earlier studies are based on wavelet analysis (WA), which can be used to recover potential information from nonstationary signals. The primary characteristic of WA is its capacity to localize a process on a temporal scale (Nourani et al. 2014). Multiple studies have shown that utilizing wavelet-decomposed time series as model inputs can produce more accurate predictions than using the original time series (Liu et al. 2011; Nourani et al. 2012; Nourani et al. 2014; Badrzadeh et al. 2016). However, several disadvantages for WA are not insignificant. For instance, poor mother wavelet functions or a low level of decomposition will reduce the accuracy of predictions.

As previously mentioned, ARIMA stands out as a popular choice for drought prediction owing to its enhanced flexibility and ability to capture comprehensive data changes over time. However, relying solely on the forecasting accuracy of a single model may fall short in meeting the demands of effective drought prediction. Moreover, precipitation data often exhibits nonlinear and nonstationary features. To enhance prediction accuracy, researchers turn to the use of hybrid models in the field of drought research.

In addition to wavelet transform, empirical mode decomposition (EMD) is one such method employed for processing nonstationary and nonlinear signals with frequency components. Unlike wavelet transformation, which may suffer from the limitations of poor mother wavelet functions or a low level of decomposition leading to lower accuracy in drought prediction, EMD emerges as a valuable technique. Combining EMD with ARIMA proves to be an effective strategy for improving time series drought prediction. This integration addresses the shortcomings associated with the purely linear and stationary nature of ARIMA, making the hybrid model better suited to handle the complexities inherent in drought forecasting (Libanda & Nkolola 2022; Xu et al. 2022).

Huang et al. (1998) proposed EMD, a data-adaptive technique that operates directly in temporal space as opposed to the equivalent frequency space as in WA (Huang & Wu 2008). Without requiring any prior basis functions, EMD sifts nonlinear and nonstationary series into several intrinsic mode functions (IMFs) and a residual component (Lee & Ouarda 2010; Zhu et al. 2016). As a data preparation technique, it can enhance the functionality of data-driven models (Wu & Huang 2009). According to Wang et al. (2018), EMD performs better than spectral time series analysis technology in analyzing the distinctive scales of annual streamflow time series. To forecast annual streamflow, Zhang et al. (2016) proposed a hybrid data-driven model that combines the external force variable, radial basis function neural networks, and EMD. They came to the conclusion that the hybrid model is workable for data-driven hydrologic forecasting in complex socio-hydrologic systems.

Wu & Huang (2009) proposed ensemble empirical mode decomposition (EEMD) with the aid of white noise to solve the mode mixing problem of the EMD method, which arises when a single IMF contains signals of widely different scales or when a signal of the same scale resides in different IMF components of the EMD method. With the addition of white noise and without the requirement for any basis functions, it can clearly combine various scale signals into suitable IMFs (Wu & Huang 2009). The mode mixing issue with the EMD approach is effectively resolved by EEMD. The genuine IMF component is defined by EEMD as the mean of an ensemble of trials, where each trial comprises the original signal and a finite-amplitude white noise. With its constituent parts belonging to various scales, the additional white noise is uniformly present over the entire time–frequency space. The mode mixing issue can be somewhat eliminated by EEMD, while the physical uniqueness of the decomposition can be maintained (Wu & Huang 2009). By employing EEMD, each component represents the real instinct change rule for a specific timescale of the original data.

For the prediction of hydrological streamflow, hybrid approaches combining decomposition and reconfiguration models have been quite popular recently. For monthly streamflow forecasting, Zhu et al. (2016) combined EMD and SVM models, and the hybrid model's best mean absolute percentage error (MAPE) is 17%. To create three hybrid models for the monthly streamflow forecasting of two stations, Zhang et al. (2015) merged three preprocessing techniques – WA, EMD, and single spectrum analysis (SSA) – with the autoregressive moving average (ARMA) model. These models are WA-ARMA, EMD-ARMA, and SSA-ARMA. EMD-ANN hybrid models with a predictive accuracy correlation coefficient (R) of 0.801 were developed by Kisi et al. (2014). For hydrological forecasting, attempts to couple EEMD with ANN were also made (Di et al. 2014; Wang et al. 2015; Zhao & Chen 2015). The original hydrological time series was broken down into several components by EEMD, and these components were used as inputs in ANN models for forecasting, and the outcomes supported its validity.

However, not many scientists have used EEMD to predict droughts. As a result, the primary goal of this study is to provide a novel hybrid method that effectively combines the benefits of two distinct models including the EEMD and ARIMA, for SPI-based drought prediction. Previous research on the drought in Afghanistan did not examine this hybrid method. The predictive hybrid model developed in this study using the SPI time series contributes to the existing knowledge and offers a fresh method for predicting future drought conditions. In summary, our research introduces a novel hybrid approach, combining EEMD with the ARIMA model for SPI-based drought prediction in Herat, Afghanistan – a region prone to severe drought but lacking sufficient scientific exploration. This study addresses a crucial literature gap by introducing EEMD, an underexplored technique in drought forecasting. By combining the strengths of ARIMA and EEMD, the hybrid model presented here aims to enhance predictive accuracy. This fills a significant gap in the literature and provides a valuable tool for water resource management in the drought-prone area of Herat, Afghanistan.

Study area and data

Afghanistan has been suffering from a severe to below-normal drought, which has impacted the accessibility and availability of clean drinking water. Around 79% of homeowners reported having insufficient access to water for basic needs like drinking, cooking, bathing, and hygiene, according to the study by Mayar (2021). Therefore, by foreseeing drought, policymakers and managers of the water supply can at least lessen its disastrous impacts by taking sensible precautions in advance. The World Bank Climate Change Knowledge Portal (https://climateknowledgeportal.worldbank.org) was used to obtain the monthly precipitation data for Herat, the largest and most significant urban area in western Afghanistan, from January 1970 to December 2019 for this study.

Afghanistan is located in the southern boundary of the Eurasian plate, which was suggested during the Permian and Triassic periods, according to Ruleman et al. (2007). It consists of three thick sedimentary rock regions: the Southeastern Katawaz region, the Northern Afghanistan Basin, and the Southwestern Afghanistan Basin. The majority of Herat City's location is within the Basin of Southwest Afghanistan. The Hari Rud-Murghab River System, a 1,100 km long fluvial system, dominates the extensively irrigated areas surrounding Herat City and accounts for about 12% of Afghanistan's total water resources. At an elevation of 922 m, Herat in western Afghanistan is located on the Hari Rud River with the latitude and longitude coordinates of 34.343044 and 62.199074, respectively. Figure 1 depicts the geographical region under investigation.

Herat province is one of the provinces in Afghanistan grappling seriously with the issue of water scarcity. Mismanagement of surface waters has led farmers to increasingly rely on groundwater resources in recent years, significantly diminishing the water tables in this province. The substantial decrease in rainfall, irrational water consumption, and unregulated digging of deep wells are major factors contributing to the fundamental decline in groundwater levels, exacerbating drought conditions in the region. Concerns among the citizens of Herat are heightened as the current water consumption solely relies on underground sources, potentially leading to severe water shortages and increased aridity in the future. Moreover, in a country like Afghanistan, characterized by arid regions, environmental threats and water resource preservation are of paramount importance, especially considering the influence reaching the northern parts of the country. The current precipitation levels, water resource limitations, and climatic conditions underscore the necessity for drought preparedness and serious measures to address and mitigate its impacts when it occurs. Drought signifies various environmental effects on living organisms, animals, and humans, whereas aridity is a generalized concept not applicable to environmental elements. Proper management and reduction of drought-related consequences in regional development require strategic planning and proactive measures, demanding the application of sufficient knowledge in drought prediction (Timuri 2023).

Standardized precipitation index

The SPI, a widely used DI for detecting drought events, was first proposed by McKee et al. (1995). It uses historical precipitation data in space and includes both positive and negative values; the former show surpluses, while the latter show shortages. Despite the fact that there are numerous DI, the WMO tried to create a standard index that can be used as a starting point for all regions and nations, as opposed to other DI that can be applied in only certain areas and require large datasets and complicated execution procedures (Yihdego et al. 2019). In contrast to other DIs such as the Palmer index and crop moisture index, the SPI is solely based on precipitation, which makes its evaluation relatively simple (Cacciamani et al. 2007). SPI is regarded as the DI for illustrating the variability in Eastern African drought (Ntale & Gan 2003), and it can describe drought at various timescales, which is one of the main advantages of using SPI to forecast the occurrence of droughts.

According to Belayneh & Adamowski (2013), SPI can be calculated by fitting a probability distribution to monthly precipitation aggregate series (3, 6, 9, 12, and 24 months). This probability density function (PDF) was then converted into a normal standardized index, whose values classify the category of drought characteristics for each location and timescale.

Gamma distribution is used to fit lengthy periods of historical precipitation records because it closely matches the sequence of precipitation episodes. The PDF derived from the gamma distribution is shown in Equation (1).
formula
(1)
where the scale, shape variables, quantity of precipitation, and gamma function are represented by , , x, and , respectively. The best estimates of and bounds are provided by Equations (2) and (3) (Guttman 1999).
formula
(2)
where the rainfall average and the total number of observations are represented by the letters x and n, respectively.
formula
(3)
formula
(4)
The variables required to calculate the cumulative probability for nonzero rainfalls are shown in Equation (5).
formula
(5)
where
For x = 0, the gamma function is uncertain, the rainfall time series data may include zero precipitation, the cumulative probability of zero and nonzero precipitation is calculated, and Equation (6) is used to determine H(x).
formula
(6)
where q is the probability of no precipitation, is the number of zeros (which occur in a rainfall time series), and is calculated using m/n.
Once the cumulative probability is changed to a standardized normal distribution (Sönmez et al. 2005), as demonstrated by Equations (7) and (8), the SPI mean and variance would equal zero and one, respectively.
formula
(7)
formula
(8)

The constants c0 = 2.515517, c1 = 0.802853, c2 = 0.010328, d1 = 1.432788, d2 = 0.189269, and d3 = 0.001308.

Table 1 displays the SPI categorization.

To sum up, the SPI was chosen for this study due to its unique advantages over other drought measures. Supported by the WMO, SPI provides robust spatial interoperability across a variety of climate zones with a simpler computing process. Because of its probabilistic properties, it is a solid option, particularly in regions with data-sparse characteristics, where immediate access to additional parameters may be challenging. Furthermore, a thorough understanding of the various types of droughts – from hydrological and socioeconomic droughts to meteorological and agricultural droughts – is made possible by SPI's versatility in computing different timescales. This study highlights SPI's critical role in tackling water resource management issues, particularly in Herat province, by utilizing its versatility to improve drought prediction in Afghanistan.

Table 1

Classification of drought and wetness using SPI

SPI valuesCategory
+2 and higher than +2 Extremely wet 
1.5 to 1.99 Very wet 
1 to 1.49 Moderately wet 
−0.99 to 0.99 Nearly normal 
−1.49 to −1 Moderate drought 
−1.99 to −1.5 Severe drought 
−2 and less than −2 Extreme drought 
SPI valuesCategory
+2 and higher than +2 Extremely wet 
1.5 to 1.99 Very wet 
1 to 1.49 Moderately wet 
−0.99 to 0.99 Nearly normal 
−1.49 to −1 Moderate drought 
−1.99 to −1.5 Severe drought 
−2 and less than −2 Extreme drought 

ARIMA model

Moving average (MA) and autoregressive (AR) models can be efficiently combined to create the broadly applicable class of time series models known as ARMA models. However, they can be utilized when the data are stationary. In the ARMA model, the present value of the time series is described as a linear aggregate of p past values and a weighted sum of q previous deviations (original value minus fitted value of previous data). By enabling differencing between data series, this class of models may be extended to nonstationary series. ARIMA models gained popularity thanks to Box and Jenkins in 1976. A model belonging to the ARIMA family is identified by three parameters (p, d, and q), which can have zero or positive integral values. The generic nonseasonal ARIMA model is AR to order p and MA to order q and works on the dth difference of the time series zt. It is possible to write the broad nonseasonal ARIMA model as follows:
formula
(9)
where and are polynomials of order p and q, respectively. The operators and will be described as follows:
formula
(10)
formula
(11)

SARIMA model

The ARIMA model was expanded by Box et al. (1994) to account for seasonality, and the result is known as the SARIMA model. By adding seasonal periodic change to a generic ARIMA mode, the SARIMA model is examined. Its denotation is ARIMA(p,d,q)(P,D,Q)S, where (p,d,q) is the model's nonseasonal component and (P,D,Q)S is its seasonal component. This may be stated in the following way:
formula
(12)
where P is the order of the seasonal AR model, Q is the order of the seasonal MA model, D is the seasonal differencing, S is the length of a season, p is the order of the nonseasonal AR model, q is the order of the nonseasonal MA model, and d is the order of the general difference.

Empirical mode decomposition

A nonstationary data processing method is EMD (Huang et al. 1998; Wang et al. 2018). Numerous applications of EMD have been made in the field of hydrological analysis to demonstrate its effectiveness without the necessity for basic functions (Lee & Ouarda 2010; Zhu et al. 2016). Different timescale features are present in hydrological time series, and these characteristics may be recovered using EMD as IMFs with various timescales. Each IMF must adhere to the following two requirements: (1) In the whole dataset, the number of zero crossings and local extreme values must be equal or separated by one. (2) At any time point, the mean values of both the upper envelope, which is defined by the local maximum, and the lower envelope, which is defined by the local lowest, must be zero.

The EMD procedure can be described as follows (Wu & Huang 2009):

  • 1.

    Let stand for the SPI original time series. To create the upper envelope, , and lower envelope, , all local extremes of x(t) are detected, and all maxima and minima are connected by a cubic spline line (Hou & Andrews 1978). The cubic spline might have a large swing at the endpoint, which is a severe issue for the spline fitting approach. To solve this issue, Huang & Wu (2008) inserted characteristic waves at the endpoints of the original time series, which are determined by the two successive extrema for the frequency and amplitude of the additional waves. Large swings have been effectively contained using this technique.

  • 2.
    The following equation illustrates how the mean of upper envelope and lower envelope is generated.
    formula
    (13)
  • 3.
    To obtain the component h(t), as stated in Equation (14), m(t) is subtracted from the initial time series x(t).
    formula
    (14)

It is verified whether that the two conditions of the IMF are met by the series h(t). If not, h(t) is replaced with x(t) and steps 1–3 are continued until h(t) satisfies the requirements of IMF, at which point h(t) becomes IMF I1(t). The residual r1(t), as shown in Equation (15), is obtained by subtracting I1(t) from the initial time series x(t).
formula
(15)
Then, steps 1–3 are repeated using r1(t) as the updating original time series to produce IMF I1(t), I2(t),…, In(t), and then the residual series rn(t), which is a monotonic function or a function with just one extreme value, and from which no more IMF can be deduced, is derived. Then, Equation (16) may be used to represent the initial time series.
formula
(16)

By using the aforementioned procedure, the initial SPI series can be broken down into a number of straightforward IMFs and a residual sequence, revealing the distinct scale and probable trend of time series. The EMD technique is an empirical, direct, and adaptable approach for data analysis that, in contrast to all the prior methods, analyzes data in temporal space rather than matching frequency space (Huang & Wu 2008).

Ensemble empirical mode decomposition

Due to signal intermittency, EMD may have mode mixing issues that cause timescale inaccuracies in IMF components, which prevents these components from accurately representing the real timescale features of the original time series (Wu & Huang 2009). The EEMD approach was put out by Wu & Huang (2009) to address such mode mixing issues. The average of a group of trials, each of which consists of a signal and a white noise with a limited amplitude, is what is referred to as the actual IMF components. The extra white noise would evenly spread the time–frequency space with a number of components at various scales. The signals with varied scales are automatically projected onto the appropriate scales of reference provided by the white noise in the background when they are added to the uniform white backdrop. Thus, this strategy effectively solves the mode mixing problem (Wu & Huang 2009). As a result, each IMF component exhibits characteristics of the original data's appropriate timescale. The following succinct summary of EEMD implementation is provided by Wu & Huang (2009):

  • 1.
    To the original series, add random white noise that follows the normal distribution, where is known and is equal to zero:
    formula
where y(t) is the initial signal, yi(t) is the signal after the ith addition of white noise, and ni(t) is the white noise.

  • 2.

    Use EMD to deconstruct time series data after adding white noise.

  • 3.
    Repeat steps 1 through 2 with various white noise, adding fresh white noise each time, integrating the results of the decomposition, and taking the average as the final result:
    formula
    (17)
    where N is indicating the number of white noises added to the model, and is the ith IMF following the addition of the jth white noise.

Because EEMD is effective at addressing the nonstationarity and nonlinearity present in drought time series data, it is frequently used as a preprocessing technique in conjunction with models such as ARIMA, ANN, SVM, and others in the field of drought forecasting. By breaking down data into IMFs, EEMD makes it easier to conduct a multiscale study and captures the different temporal scales that impact drought events. Its versatility in handling intricate patterns and capacity to boost the signal-to-noise ratio render it advantageous in distinguishing meaningful signals from background interference, hence augmenting the precision of ensuing forecasting models. Enhancing forecasting performance and supporting decision-making in water resource management and related applications, EEMD offers a more nuanced depiction of the underlying dynamics and mitigates assumptions violations in standard models (Başakın et al. 2021; Liu et al. 2022).

By decomposing the original data series, EEMD can effectively handle nonlinear and nonstationary data series and extract valuable information from the data. Wang et al. (2015) present the use of EEMD in conjunction with the ARIMA model to forecast yearly runoff time series. The findings suggest that EEMD can significantly improve predicting accuracy. Because of its strong generalization capabilities, the radial basis function neural network (RBFN) model is frequently employed in the time series analysis. It is capable of effectively forecasting data series. In recent years, ANN has been applied extensively and with good results in a wide range of fields. To solve the hydrologic problem, a hybrid model called EEMD–RBFN–SVM was presented. It was based on EEMD, RBFN, and SVM. The EEMD–RBFN–SVM model's performance was assessed using precipitation data from the Qinghai Yushu Tibetan Autonomous Region, which exhibits clear outperformance of hybrid EEMD with other RBFN and SVM in the context of hydrological time series prediction (Jiao et al. 2016).

The success of combining EEMD with various statistical and machine learning methods to tackle nonlinear and nonstationary time series has strongly encouraged us to incorporate this preprocessing approach into the traditional ARIMA framework for drought forecasting.

2.7. Hybrid EEMD-ARIMA model

Since the hybrid model combines multiple approaches, it is expected to perform better than other SPI-based drought prediction models. The hybrid method utilizes the complementing strengths of various models by combining preprocessing techniques such as EEMD with statistical methods like ARIMA. To address the nonlinear and nonstationary properties included in SPI time series data, this integration is especially important because it enables a more thorough description of intricate patterns. Furthermore, the hybrid model is expected to outperform in capturing and forecasting complex dynamics of drought conditions based on SPI metrics due to its adaptability to varied temporal scales, decreased sensitivity to model assumptions, and improved accuracy through the collaborative expertise of different methods (Wang et al. 2018; Alawsi et al. 2022).

In this study, the original SPI time series is first divided into multiple components with distinct timescales by EEMD to obtain stationary subsequences. Each of these subsequences is then predicted one time step ahead by ARIMA, and the sum of these predicted subsequences is taken as the prediction value of the original SPI for the following time step. The hybrid EEMD-ARIMA models' flowchart is shown in Figure 2, and the stages for implementation are as follows:
  • 1.

    First, using the EEMD methodology, the original yearly SPI time series y(t) (t = 1, 2,…, n) is divided into m IMF components, ci(t), i = 1, 2,…, m, and one residual component rm(t).

  • 2.

    Second, the ARIMA model is applied as a forecasting tool to model the decomposed IMF and residual components and produce the corresponding prediction for each extracted IMF component and residual component, respectively.

  • 3.

    Finally, the forecasting outcomes of all extracted IMF and residual components acquired by ARIMA are aggregated to provide an output that may be utilized as the final forecasting outcome for the initial yearly SPI time series.

Figure 2

EEMD-ARIMA framework.

Figure 2

EEMD-ARIMA framework.

Close modal

In conclusion, the hybrid EEMD-ARIMA forecasting model implements the decomposition and ensemble concepts. The ensemble creates a consensus forecast based on the original data, whereas the decomposition simplifies the forecasting work. Herat annual precipitation time series from January 1970 to December 2019 are used for testing purposes in this study to verify and make the pattern of the extracted IMFs and residual components reflect the forecasting model and to improve forecasting performance.

Evaluation indicators

By contrasting observation and prediction, different models' predictive abilities are measured. The performance of hybrid models, as stated in Equations (18)–(21), is estimated using the RMSE, MAE, and MAPE. While MAE calculates the fitness of moderate streamflow from a more well-rounded viewpoint, RMSE may determine the degree of fitting between the anticipated and observed data with a high value and MAPE is used to assess the fitness between the predicted and the observed data with a moderate value (Zhu et al. 2016). R2 indicates the proportion of variability in observed data accounted for by the model. Further information about these performance evaluation metrics are available in the literature (Moriasi et al. 2007; Tongal 2013; Buyukyildiz et al. 2014; Wu et al. 2021; Koycegiz & Buyukyildiz 2019, 2022, 2023).
formula
(18)
formula
(19)
formula
(20)
formula
(21)
where n is the number of data, Yi indicates the observed data, shows the mean of observed data, and represents the predicted data. Based on these criteria, a better model performance will be gained with smaller RMSE and MAE, and MAPE and large R2.

Key metrics such RMSE, MAE, MAPE, and R2 have specific meanings in the context of drought prediction and management. Decision-makers can better comprehend and gauge the severity of the situation when the forecasting model more correctly captures the magnitude of drought occurrences, as indicated by the reduced RMSE and MAE. Simultaneously, a lower MAPE signifies a closer percentage alignment between the expected and observed values, providing reliable data for mitigation measures and resource allocation. Furthermore, a higher R2 indicates that a greater percentage of the variability in the observed drought conditions can be explained by the model, providing confidence in its capacity to replicate drought patterns and enabling sound resource allocation and management decision-making.

R2 is a significant parameter to measure the explanatory power of the model and its capacity to explain the variance in the observed SPI values in the context of SPI-based drought forecasting. A higher R2 indicates that the forecasting model can account for a greater share of the variability in SPI, which represents standardized precipitation anomalies. Practically speaking, a high R2 indicates that the model is good at identifying and recreating the dynamics and patterns present in SPI time series data. This reflects the degree to which the model offers a trustworthy depiction of observed drought conditions, which is important information for decision-makers in drought management. Conversely, a low R2 can suggest that the model is not doing a good enough job of explaining the variability in SPI, which could result in fewer accurate forecasts. Consequently, a careful analysis of R2 in SPI-based drought forecasting highlights the explanatory ability of the model, affecting the level of trust that managers of water resources might have in the forecasted results and facilitating better-informed decision-making in drought-prone areas.

In this study, the objective is to enhance the accuracy of the ARIMA model in predicting drought events in Herat province, Afghanistan, using the EEMD-ARIMA model. The dataset is divided into two periods: a training phase from 1970 to 2009 and a validation phase from 2010 to 2019. The model development process consists of three main stages: model identification, estimation of unknown parameters, and diagnostic checking. The original SPI values are decomposed using EEMD, and R software is employed for model development and future drought prediction. To categorize drought events based on their duration, different SPI values are calculated using total precipitation data for running periods of 3, 6, 9, and 12. Specifically, SPI 3 indicates short-term drought, while SPI 6 and 9 represent mid-term drought events. SPI 12 is utilized to characterize long-term drought. Figure 3 depicts the plotted values of different SPIs based on historical records. The figure demonstrates a noticeable oscillation within each SPI.
Figure 3

SPI 3, 6, 9, and 12 plots based on historical data from Herat (the case study area).

Figure 3

SPI 3, 6, 9, and 12 plots based on historical data from Herat (the case study area).

Close modal

Decomposing SPI time series using EEMD

White noise magnitude and ensemble size are two critical parameters that directly affect the performance of decomposing data, and they should be fixed before employing EEMD procedure. The well-established statistical rule for fixing these two parameters was described by Wu & Huang (2009), which is defined as follows:
formula
(22)

where n represents ensemble number, is the amplitude of added white noise, and is indicator of the standard deviation of error. Wu & Huang (2009) came to the conclusion that the effect of the additional white noise should diminish using Equation (21) and noted that if the additional noise amplitude is too small, it might not result in the change of extrema that the EMD depends on; in contrast, if the additional noise is too large, it would produce redundant IMF components. According to the experimental findings in the majority of cases, the injected noise's magnitude should be roughly 0.2 times the standard deviation of the sample data (Wu & Huang 2009). Zhang et al. (2010) studied the EEMD parameter settings for various applications, and their results on the additional white noise amplitude are consistent with the results presented by Wu & Huang (2009). The EEMD approach is used in this work with an ensemble size of 100 and a white noise amplitude of 0.2 times standard deviation. Wu & Huang (2009) provide a detailed examination of the impacts of noise and ensemble averaging on decomposition with reference to these factors; thus, it is not necessary to repeat it here.

Hence, this study follows the statistical principles outlined by Wu & Huang (2009) to determine the amplitude of added white noise, ensemble number, and the indicator of the standard deviation of error (σ). These principles guide the selection of white noise magnitude and ensemble size for the EEMD process. The parameters are fixed based on Equation (22). Aiming for an ideal balance, the selected white noise amplitude was roughly 0.2 times the sample data standard deviation. This allowed for enough impact without adding unnecessary information, which is consistent with the results of the experiments and the literature (Zhang et al. 2010). The ensemble size, set at 100, followed established practices, considering the robustness of the decomposition process for precipitation data. So, to ensure the physical relevance and accuracy of the EEMD decomposition in the context of drought forecasting for SPI time series, established guidelines and literature findings were taken into consideration while choosing the values.

As stated earlier, each SPI time series can be decomposed into a number of separate IMFs and a single residue using the EEMD technique. Each original SPI (SPI 3, 6, 9, and 12) time series are divided into eight distinct IMF components in the order of decreasing frequency and one residue component, respectively. The outcomes for SPI 12 are depicted in Figure 4 as an example.
Figure 4

Decomposition of SPI 12 time series using EEMD.

Figure 4

Decomposition of SPI 12 time series using EEMD.

Close modal

All IMFs were found to have unique amplitudes, wavelengths, and frequencies. The first IMF had the smallest wavelength, highest frequency, and greatest amplitude. In comparison to the first IMF, the second IMF had a lesser amplitude, lower frequency, and longer wavelengths. Similar to this, the pattern continues to the last IMF in which the wavelength grows while the matching frequency and amplitude decreases. As a result, the decomposition offered by EEMD is physically relevant with various scales for each independent IMF (Debert et al. 2011). The residue component of the time series for monthly precipitation shows the overall trend. Decomposition can be therefore used to convert nonstationary, nonlinear time series into stationary time series and to enhance drought forecasting accuracy.

ARIMA modeling approach

In general, the ARIMA model is frequently used for time series forecasting for the Box–Jenkins approach (Box & Jenkins 1976). However, for the ARIMA model to be applied, the time series must be stationary. In other words, the ARIMA algorithm assumes that the process maintains equilibrium at a fixed mean level. So, checking the stationarity of the data is the first step in time series approaches. Plots of the autocorrelation function (ACF) and partial autocorrelation function (PACF) can be used to check data stationarity, which a rapid fluctuation in time series data indicate nonstationarity. The nonstationary nature of the time series is shown by a slow decaying ACF plot, and it is eliminated using differencing treatment to produce stationary data (Box et al. 1994). In this study, to check the stationarity of monthly original SPI and decomposed SPI series augmented Dickey–Fuller (ADF) is utilized (Elliott et al. 1996).

Model identification, parameter estimation, and diagnosis checking are the three primary phases that make up the ARIMA model, as we are all aware. The identification step, which consists of two stages, is crucial among these three procedures: (1) if necessary, proper series differencing is carried out to attain stationarity and normality and (2) the order of the AR and MA parts of ARIMA model is identified. The main tools used by Box and Jenkins (Box & Jenkins 1976) to determine the order of the ARIMA model were the ACF and PACF of the sample data. If sample data is an AR(p) model, the PACF cutoff is at lag p. In contrast, if sample data are an MA(q) process, the ACF has a cutoff at lag q. However, when dealing with mixed ARIMA processes, the PACF and ACF technique is not very helpful. In most cases, a quick glance at the ACF and PACF graphs would not reveal the values of p and q for mixed models (Chan 1999). The maximum likelihood approach, which was put forth by Box & Jenkins (1976), was used in this study to estimate the ARIMA model's parameters, which is a technique for determining a statistical model's parameters.

The Akaike's information criterion (AIC) (Shibata 1976), the Bayes information criterion (BIC) (Schwarz 1978), the final prediction error criterion (Akaike 1969), and other information-theoretic approach-based identification techniques have been presented. The best fitted model is chosen in this study based on AIC and BIC. The model's parameters must be estimated once an acceptable model has been selected. The mathematical formulation of these criteria is described as follows:
formula
(23)
formula
(24)
where k is illustrating the number of parameters in the model , L is the likelihood function of the ARIMA model, and n is the number of data in the model. The goal of model identification is to uncover a fitting subclass of the ARIMA model that accurately fits the time series data.

Proposed EEMD-ARIMA approach

According to the aforementioned description, the ADF test is first applied to the original monthly SPI time series and all deconstructed SPI components. Table 2 provides a summary of the outcomes of the ADF unit root testing. As seen in Table 2, for SPI 3, both the original and decomposed IMF 1–IMF 5 time series are stationary. The ADF test indicated that IMF 8 and the residue became stationary after applying one differencing. In addition, IMF 6 and IMF 7 achieved stationarity after applying a differencing order of 2. In the case of SPI 6, the original and decomposed IMF 1–IMF 5 are already stationary. In contrast, IMF 6 became stationary by applying first-order differencing, while IMF 7 and IMF 8 achieved stationarity through the application of second-order differencing. Original IMF 1–IMF 4 and residue of SPI 9 are inherently stationary. To achieve stationary series for IMF 5–IMF 8, the second-order differencing is employed. The original data, along with IMF 1 to IMF 4 and the residue of SPI 9, is inherently stationary. To achieve stationary series for IMF 5–IMF 8, second-order differencing is applied. In a similar manner, for SPI 12, to attain stationary series for IMF 5–IMF 8 as well as the residue, second-order differencing is applied. However, the ADF test for the original data and IMF 1–IMF 8 has already confirmed the presence of stationary series, negating the need for applying differencing.

Table 2

Augmented Dickey–Fuller test for stationarity in SPI series

Augmented Dickey–Fuller test for SPI series
SPIi-Statisticp-ValueSPIt-Statisticp-Value
SPI 3 Original −8.7072 0.01 SPI 6 Original −6.7912 0.01 
IMF 1 −11.312 0.01 IMF 1 −11.044 0.01 
IMF 2 −11.109 0.01 IMF 2 −10.817 0.01 
IMF 3 −10.568 0.01 IMF 3 −11.866 0.01 
IMF 4 −6.7024 0.01 IMF 4 −4.7236 0.01 
IMF 5 −5.0946 0.01 IMF 5 −5.2543 0.01 
IMF 6 −4.2457 0.347 IMF 6 −3.7869 0.0197 
IMF 7 −3.1274 0.1011 IMF 7 −2.1965 0.4952 
IMF 8 2.5335 0.99 IMF 8 1.4107 0.99 
Residue 0.22081 0.99 Residue −193117 0.01 
SPI 9 Original −8.2453 0.01 SPI 12 Original −6.2147 0.01 
IMF 1 −10.449 0.01 IMF 1 −10.319 0.01 
IMF 2 −11.989 0.01 IMF 2 −13.536 0.01 
IMF 3 −9.9568 0.01 IMF 3 −12.449 0.01 
IMF 4 −7.0845 0.01 IMF 4 −7.4394 0.01 
IMF 5 −1.6962 0.7069 IMF 5 −1.4495 0.8113 
IMF 6 −3.8991 0.0141 IMF 6 −3.7497 0.0215 
IMF 7 −2.0814 0.5439 IMF 7 −1.7361 0.69 
IMF 8 1.3193 0.99 IMF 8 2.6928 0.99 
Residue 0.54257 0.99 Residue 0.067147 0.99 
Augmented Dickey–Fuller test for SPI series
SPIi-Statisticp-ValueSPIt-Statisticp-Value
SPI 3 Original −8.7072 0.01 SPI 6 Original −6.7912 0.01 
IMF 1 −11.312 0.01 IMF 1 −11.044 0.01 
IMF 2 −11.109 0.01 IMF 2 −10.817 0.01 
IMF 3 −10.568 0.01 IMF 3 −11.866 0.01 
IMF 4 −6.7024 0.01 IMF 4 −4.7236 0.01 
IMF 5 −5.0946 0.01 IMF 5 −5.2543 0.01 
IMF 6 −4.2457 0.347 IMF 6 −3.7869 0.0197 
IMF 7 −3.1274 0.1011 IMF 7 −2.1965 0.4952 
IMF 8 2.5335 0.99 IMF 8 1.4107 0.99 
Residue 0.22081 0.99 Residue −193117 0.01 
SPI 9 Original −8.2453 0.01 SPI 12 Original −6.2147 0.01 
IMF 1 −10.449 0.01 IMF 1 −10.319 0.01 
IMF 2 −11.989 0.01 IMF 2 −13.536 0.01 
IMF 3 −9.9568 0.01 IMF 3 −12.449 0.01 
IMF 4 −7.0845 0.01 IMF 4 −7.4394 0.01 
IMF 5 −1.6962 0.7069 IMF 5 −1.4495 0.8113 
IMF 6 −3.8991 0.0141 IMF 6 −3.7497 0.0215 
IMF 7 −2.0814 0.5439 IMF 7 −1.7361 0.69 
IMF 8 1.3193 0.99 IMF 8 2.6928 0.99 
Residue 0.54257 0.99 Residue 0.067147 0.99 

As a supportive visualization method, we also used ACF and partial autocorrelation plots to check the stationarity of SPI series. To determine stationarity, a rapid decay is observed in ACF, a sudden cutoff takes place in PACF, and it is ensured that there are no persistent spikes in either plot beyond specific lags. On the other hand, if there is a slow decay in ACF or persistent spikes in PACF, it might indicate nonstationarity, meaning the data have a trend or repeating pattern. Overall, we are checking for patterns that show the data do not have a persistent trend for a stationary series. Combining both ACF and PACF helps us get a clear picture.

In this study, Figures 5 and 6 depict the ACF and PACF for SPI3 and selected decomposed components (IMF1, IMF6, and IMF8) as illustrative examples. Figure 5 reveals that the original SPI3 and decomposed IMF1 exhibit stationarity. Conversely, in Figure 6, the ACF and PACF plots for IMF6 and IMF8 suggest nonstationary series. These findings align with the results obtained from the ADF test.
Figure 5

ACF and PACF plot of SPI 3 and decomposed component IMF1.

Figure 5

ACF and PACF plot of SPI 3 and decomposed component IMF1.

Close modal
Figure 6

ACF and PACF plots of decomposed components IMF6 and IMF8 in SPI 3.

Figure 6

ACF and PACF plots of decomposed components IMF6 and IMF8 in SPI 3.

Close modal

To determine the order of the ARIMA model, we utilized both ACF and PACF plots, in addition to the ‘auto.arima’ function in R. First, it is essential to confirm the stationarity of the data using the ADF test. If the series is found to be nonstationary, differencing is applied to achieve stationarity. Next, determine the order of auto regressive (p) and MA (q) terms by looking at the ACF and PACF plots of the differenced series. The PACF plot, which shows partial autocorrelations, assists in detecting p order, while the ACF plot, which shows autocorrelations at different lags, assists in recognizing the necessity for q terms. Examine both plots for notable spikes and cutoff patterns. A severe cutoff or notable spike at lag k in the PACF figure indicates that a p term of order k is required. Likewise, an order of q is shown if the ACF plot stops at lag k. To prevent overdifferencing and to look for unit roots, care should be given and model orders adjusted appropriately. The aim is to achieve a balance between model complexity and efficacy in identifying the underlying patterns in the time series.

The best-selected model for each original SPI and its decomposed components was chosen based on the AIC and BIC criteria. The models with the minimum AIC and BIC were considered the best choices among the other candidates. Table 3 indicates the best-selected ARIMA/SARIMA models for each SPI and their corresponding decomposed counterparts, while Table 4 provides the estimated parameters for these selected models.

Table 3

Summary of the best-selected ARIMA/SARIMA model for each decomposed and original SPI se based on AIC and BIC criterion

SPI seriesModelAICBIC
SPI 3 Original ARIMA(2,0,2) 1168.39 1193.41 
IMF 1 ARIMA(2,0,1)(2,0,1)[3] 583.16 612.37 
IMF 2 ARIMA(3,0,3)(1,0,1)[3] −683.83 −642.09 
IMF 3 ARIMA(4,0,3) −2911.12 −2873.56 
IMF 4 ARIMA(3,0,4)(1,0,1)[3] −5900.84 −5854.93 
IMF 5 ARIMA(5,0,4)(1,0,1)[3] −9230.57 −9176.31 
IMF 6 ARIMA(1,2,3)(0,0,1)[3] −10,817.71 −10,792.69 
IMF 7 ARIMA(2,2,1) −13,175.09 −13,158.41 
IMF 8 ARIMA(1,2,3) −16,675.75 −16,654.9 
Residue ARIMA(2,2,2) −21,672.52 −21,651.67 
SPI 6 Original ARIMA(2,0,1)(2,0,0)[6] 1035.7 1064.85 
IMF 1 ARIMA(2,0,2)(2,0,0)[6] 401.46 430.68 
IMF 2 ARIMA(2,0,4) −934.82 −901.43 
IMF 3 ARIMA(4,0,1)(2,0,2)[6] −2929.65 −2887.91 
IMF 4 ARIMA(4,0,1)(2,0,1)[6] −5443.85 −5402.11 
IMF 5 ARIMA(2,0,2)(2,0,1)[6] 404.55 442.12 
IMF 6 ARIMA(2,1,2) −10,366.43 −10,345.57 
IMF 7 ARIMA(1,2,3) −12,354.06 −12,333.21 
IMF 8 ARIMA(1,2,3) −15,755.59 −15,734.74 
Residue ARIMA(0,2,4)(0,0,2)[6] −14,086.88 −14,057.69 
SPI 9 Original ARIMA(1,0,0)(2,0,1)[9] 691.37 716.32 
IMF 1 ARIMA(6,0,6)(1,0,1)[9] −120.37 −53.59 
IMF 2 ARIMA(5,0,5)(1,0,0)[9] −1328.13 −1273.87 
IMF 3 ARIMA(5,0,3) −3384.73 −3342.99 
IMF 4 ARIMA(4,0,2)(0,0,1)[9] −5956.57 −5919 
IMF 5 ARIMA(3,2,1)(0,0,1)[9] −9195.8 −9170.78 
IMF 6 ARIMA(1,2,1)(0,0,1)[9] −9976.5 −9959.82 
IMF 7 ARIMA(1,2,2) −12,658.72 −12,642.04 
IMF 8 ARIMA(1,2,2) −15,963.44 −15,946.76 
Residue ARIMA(1,2,0) −14,610.08 −14,601.74 
SPI 12 Original ARIMA(1,1,0)(1,0,1)[12] 238.8 255.39 
IMF 1 ARIMA(4,0,5)(1,0,2)[12] −636.18 −577.74 
IMF 2 ARIMA(5,0,5)(0,0,1)[12] −1844.78 −1790.52 
IMF 3 ARIMA(6,0,6)(1,0,1)[12] −3868.3 −3801.52 
IMF 4 ARIMA(4,0,5)(1,0,3)[12] −6301.14 −6238.53 
IMF 5 ARIMA(2,2,3) −9408.67 −9383.65 
IMF 6 ARIMA(2,2,2)(0,0,1)[12] −10,505.03 −10,480.01 
IMF 7 ARIMA(1,2,2) −12,844.36 −12,827.68 
IMF 8 ARIMA(1,2,3) −16,132.19 −16,111.34 
Residue ARIMA(1,2,0) −32,363.44 −32,355.1 
SPI seriesModelAICBIC
SPI 3 Original ARIMA(2,0,2) 1168.39 1193.41 
IMF 1 ARIMA(2,0,1)(2,0,1)[3] 583.16 612.37 
IMF 2 ARIMA(3,0,3)(1,0,1)[3] −683.83 −642.09 
IMF 3 ARIMA(4,0,3) −2911.12 −2873.56 
IMF 4 ARIMA(3,0,4)(1,0,1)[3] −5900.84 −5854.93 
IMF 5 ARIMA(5,0,4)(1,0,1)[3] −9230.57 −9176.31 
IMF 6 ARIMA(1,2,3)(0,0,1)[3] −10,817.71 −10,792.69 
IMF 7 ARIMA(2,2,1) −13,175.09 −13,158.41 
IMF 8 ARIMA(1,2,3) −16,675.75 −16,654.9 
Residue ARIMA(2,2,2) −21,672.52 −21,651.67 
SPI 6 Original ARIMA(2,0,1)(2,0,0)[6] 1035.7 1064.85 
IMF 1 ARIMA(2,0,2)(2,0,0)[6] 401.46 430.68 
IMF 2 ARIMA(2,0,4) −934.82 −901.43 
IMF 3 ARIMA(4,0,1)(2,0,2)[6] −2929.65 −2887.91 
IMF 4 ARIMA(4,0,1)(2,0,1)[6] −5443.85 −5402.11 
IMF 5 ARIMA(2,0,2)(2,0,1)[6] 404.55 442.12 
IMF 6 ARIMA(2,1,2) −10,366.43 −10,345.57 
IMF 7 ARIMA(1,2,3) −12,354.06 −12,333.21 
IMF 8 ARIMA(1,2,3) −15,755.59 −15,734.74 
Residue ARIMA(0,2,4)(0,0,2)[6] −14,086.88 −14,057.69 
SPI 9 Original ARIMA(1,0,0)(2,0,1)[9] 691.37 716.32 
IMF 1 ARIMA(6,0,6)(1,0,1)[9] −120.37 −53.59 
IMF 2 ARIMA(5,0,5)(1,0,0)[9] −1328.13 −1273.87 
IMF 3 ARIMA(5,0,3) −3384.73 −3342.99 
IMF 4 ARIMA(4,0,2)(0,0,1)[9] −5956.57 −5919 
IMF 5 ARIMA(3,2,1)(0,0,1)[9] −9195.8 −9170.78 
IMF 6 ARIMA(1,2,1)(0,0,1)[9] −9976.5 −9959.82 
IMF 7 ARIMA(1,2,2) −12,658.72 −12,642.04 
IMF 8 ARIMA(1,2,2) −15,963.44 −15,946.76 
Residue ARIMA(1,2,0) −14,610.08 −14,601.74 
SPI 12 Original ARIMA(1,1,0)(1,0,1)[12] 238.8 255.39 
IMF 1 ARIMA(4,0,5)(1,0,2)[12] −636.18 −577.74 
IMF 2 ARIMA(5,0,5)(0,0,1)[12] −1844.78 −1790.52 
IMF 3 ARIMA(6,0,6)(1,0,1)[12] −3868.3 −3801.52 
IMF 4 ARIMA(4,0,5)(1,0,3)[12] −6301.14 −6238.53 
IMF 5 ARIMA(2,2,3) −9408.67 −9383.65 
IMF 6 ARIMA(2,2,2)(0,0,1)[12] −10,505.03 −10,480.01 
IMF 7 ARIMA(1,2,2) −12,844.36 −12,827.68 
IMF 8 ARIMA(1,2,3) −16,132.19 −16,111.34 
Residue ARIMA(1,2,0) −32,363.44 −32,355.1 
Table 4

Extracted structure and parameters of ARIMA/SARIMA based on AIC and BIC

Structure
SPI Intercept
SPI 3 Original −0.2409, 0.3044 0.8790, 0.2809 – – −0.0218 
IMF 1 −0.4174, −0.3940 0.0660 0.01692, 0.2351 −0.6360 – 
IMF 2 1.7332, −1.2695, 0.3661 0.5420, −0.7911, −0.5352 0.9352 −0.9888 −0.0057 
IMF 3 3.4687, −4.7641, 3.0574, −0.7772 0.0405, 0.5438, −0.4197 – – 0.0063 
IMF 4 2.6841, −2.4324, 0.7386 1.8703, 1.9749, 1.1563, 0.3704 – – −0.0022 
IMF 5 2.8356, −1.8005, −1.6866, 2.4277, −0.7766 2.2394, 2.1006, 1.0899, 0.2477 0.2829 −0.0640 −0.0080 
IMF 6 0.9511 2.1356, 2.0102, 0.7603, −0.7928 – 0.7963 – 
IMF7 1.9441, −0.9454 −0.7927 – – – 
IMF 8 0.9935 1.8167, 1.5234, 0.5849 – – – 
Residue 0.0034, 0.0997 0.0997, 0.0088 – – – 
SPI 6 Original −0.1528, 0.6894 0.9397 −0.1600, −0.1505 – −0.0081 
IMF 1 0.6380, −0.5642 −0.9607, 0.5261 −0.2298, −0.1888 – – 
IMF 2 1.1709, −0.5340 1.0466, −0.2284, −0.7986, −0.3978 – – 0.0031 
IMF 3 3.2550, −4.3106, 2.7282, −0.7098 0.5090 −0.3478, −0.5471 0.4456, 0.6053 – 
IMF 4 3.6842, −5.2066, 3.3439, −0.8250 0.6847 0.8261, −0.1560 −0.6867 −0.0064 
IMF 5 0.6725, −0.5757 −0.9951, 0.5513 −0.0891, −0.1643 −0.1556 0.0046 
IMF 6 1.9941, −0.9958 1.7473, 0.9697 – – – 
IMF7 0.9739 1.7212, 1.4532, 0.5938 – – – 
IMF 8 0.9937 2.2815, 2.0890, 0.7687 – – – 
Residue – 3.0428, 4.0657, 2.9476, 0.9383 – 1.1962, 0.8396 – 
SPI 9 Original 0.9140 – 0.3717, 0.1644 −0.7952 0.0126 
IMF 1 0.0733, −0.2958, 0.3445, 0.0454, 0.7514, −0.3630 −0.5105, 0.1553, −0.3859, −0.0024, −0.6785, 0.7361 0.1612 −0.4871 0.0017 
IMF 2 1.4464, −0.5645, −0.6913, 0.8231, −0.3078 0.4965, −0.9922, −0.7617, 0.14756, 0.3077 −0.1392 – 0.0009 
IMF 3 4.2670, −7.5117, 6.8128, −3.1820, 0.6116 −0.2681, −0.5230, −0.1493 – – 0.0026 
IMF 4 3.7167, −5.2682, 3.3759, −0.8263 1.0392, 0.3307 – 0.0520 −0.0146 
IMF 5 2.4858, −2.0056, 0.5149 0.6004 −0.0069 – – 
IMF 6 0.9895 0.9747 – 0.7081 – 
IMF 7 0.9854 1.2452, 0.6115 – – – 
IMF 8 0.9957 1.9015, 0.9485 – – – 
Residue 0.9882 – – – – 
SPI 12 Original 0.0885 – −0.2455 −0.5738 – 
IMF 1 −0.3121, 0.2519, −0.3860, −0.7530 −0.2130, −0.5854, 0.6124, 0.7488, −0.5729 −0.0138 −0.5002, 0.2241 0.0051 
IMF 2 2.5759, −3.0621, 1.9887, −0.7371, 0.1102 −0.3409, −0.8279, −0.0327, 0.3896, 0.1478 – −0.5601 0.0029 
IMF 3 1.3862, 0.7032, −1.9933, 0.3425, 1.0202, −0.5654 2.6446, 3.1023, 2.6853, 2.2607, 1.4044, 0.4271 −0.0372 −0.3341 0.0061 
IMF 4 3.7954, −5.4732, 3.5543, −0.8776 1.1042, 0.3522, −0.1246, −0.2226, −0.1185 −0.9165 0.6506, −0.0930, 0.1381 0.0004 
IMF 5 1.9686, −0.9774 1.2736, 0.6994, 0.2105 – – – 
IMF 6 1.9718, −0.9735 0.1499, −0.6899 – −0.5643 – 
IMF7 0.9880 1.7448, 0.8727 – – – 
IMF 8 0.9933 2.0196, 1.7685, 0.6679 – – – 
Residue 1.0000 – – – – 
Structure
SPI Intercept
SPI 3 Original −0.2409, 0.3044 0.8790, 0.2809 – – −0.0218 
IMF 1 −0.4174, −0.3940 0.0660 0.01692, 0.2351 −0.6360 – 
IMF 2 1.7332, −1.2695, 0.3661 0.5420, −0.7911, −0.5352 0.9352 −0.9888 −0.0057 
IMF 3 3.4687, −4.7641, 3.0574, −0.7772 0.0405, 0.5438, −0.4197 – – 0.0063 
IMF 4 2.6841, −2.4324, 0.7386 1.8703, 1.9749, 1.1563, 0.3704 – – −0.0022 
IMF 5 2.8356, −1.8005, −1.6866, 2.4277, −0.7766 2.2394, 2.1006, 1.0899, 0.2477 0.2829 −0.0640 −0.0080 
IMF 6 0.9511 2.1356, 2.0102, 0.7603, −0.7928 – 0.7963 – 
IMF7 1.9441, −0.9454 −0.7927 – – – 
IMF 8 0.9935 1.8167, 1.5234, 0.5849 – – – 
Residue 0.0034, 0.0997 0.0997, 0.0088 – – – 
SPI 6 Original −0.1528, 0.6894 0.9397 −0.1600, −0.1505 – −0.0081 
IMF 1 0.6380, −0.5642 −0.9607, 0.5261 −0.2298, −0.1888 – – 
IMF 2 1.1709, −0.5340 1.0466, −0.2284, −0.7986, −0.3978 – – 0.0031 
IMF 3 3.2550, −4.3106, 2.7282, −0.7098 0.5090 −0.3478, −0.5471 0.4456, 0.6053 – 
IMF 4 3.6842, −5.2066, 3.3439, −0.8250 0.6847 0.8261, −0.1560 −0.6867 −0.0064 
IMF 5 0.6725, −0.5757 −0.9951, 0.5513 −0.0891, −0.1643 −0.1556 0.0046 
IMF 6 1.9941, −0.9958 1.7473, 0.9697 – – – 
IMF7 0.9739 1.7212, 1.4532, 0.5938 – – – 
IMF 8 0.9937 2.2815, 2.0890, 0.7687 – – – 
Residue – 3.0428, 4.0657, 2.9476, 0.9383 – 1.1962, 0.8396 – 
SPI 9 Original 0.9140 – 0.3717, 0.1644 −0.7952 0.0126 
IMF 1 0.0733, −0.2958, 0.3445, 0.0454, 0.7514, −0.3630 −0.5105, 0.1553, −0.3859, −0.0024, −0.6785, 0.7361 0.1612 −0.4871 0.0017 
IMF 2 1.4464, −0.5645, −0.6913, 0.8231, −0.3078 0.4965, −0.9922, −0.7617, 0.14756, 0.3077 −0.1392 – 0.0009 
IMF 3 4.2670, −7.5117, 6.8128, −3.1820, 0.6116 −0.2681, −0.5230, −0.1493 – – 0.0026 
IMF 4 3.7167, −5.2682, 3.3759, −0.8263 1.0392, 0.3307 – 0.0520 −0.0146 
IMF 5 2.4858, −2.0056, 0.5149 0.6004 −0.0069 – – 
IMF 6 0.9895 0.9747 – 0.7081 – 
IMF 7 0.9854 1.2452, 0.6115 – – – 
IMF 8 0.9957 1.9015, 0.9485 – – – 
Residue 0.9882 – – – – 
SPI 12 Original 0.0885 – −0.2455 −0.5738 – 
IMF 1 −0.3121, 0.2519, −0.3860, −0.7530 −0.2130, −0.5854, 0.6124, 0.7488, −0.5729 −0.0138 −0.5002, 0.2241 0.0051 
IMF 2 2.5759, −3.0621, 1.9887, −0.7371, 0.1102 −0.3409, −0.8279, −0.0327, 0.3896, 0.1478 – −0.5601 0.0029 
IMF 3 1.3862, 0.7032, −1.9933, 0.3425, 1.0202, −0.5654 2.6446, 3.1023, 2.6853, 2.2607, 1.4044, 0.4271 −0.0372 −0.3341 0.0061 
IMF 4 3.7954, −5.4732, 3.5543, −0.8776 1.1042, 0.3522, −0.1246, −0.2226, −0.1185 −0.9165 0.6506, −0.0930, 0.1381 0.0004 
IMF 5 1.9686, −0.9774 1.2736, 0.6994, 0.2105 – – – 
IMF 6 1.9718, −0.9735 0.1499, −0.6899 – −0.5643 – 
IMF7 0.9880 1.7448, 0.8727 – – – 
IMF 8 0.9933 2.0196, 1.7685, 0.6679 – – – 
Residue 1.0000 – – – – 

Once the parameters of the selected models are estimated, we must to examine the residuals of the model to verify that the models are adequate for time series. Using a correlogram, residuals' independence was examined. Figure 7 illustrates the ACF and PACF of the original SPI 3 and SPI 12 as well as the residual ACF (RACF) and residual PACF (RPACF). The findings indicate that the majority of the RACF and RPACF values were within the 95% confidence intervals, suggesting that there was no significant association between the residuals. The confidence intervals in this context pertain to statistical significance, reinforcing the conclusion that the residuals lack a substantial correlation.
Figure 7

The ACF, RACF, PACF, and RPACF of SPI 3 and SPI 12.

Figure 7

The ACF, RACF, PACF, and RPACF of SPI 3 and SPI 12.

Close modal

The performances of ARIMA models based on the original SPI time series and the proposed EEMD-ARIMA models are assessed using the four quantitative standard statistical performance evaluation measures RMSE, MAE, MAPE, and R2. Table 5 summarizes the statistical outcomes of the various models in both training and validation periods. The most accurate models are those with the lowest RMSE, MAE, and MAPE and the highest R2.

Table 5

Forecasting evaluation indices of all SPI series

NameModelTraining
Testing
RMSEMAEMAPER2RMSEMAEMAPER2
SPI 3 ARIMA 0.811 0.632 271.170 0.380 0.658 0.510 387.054 0.353 
EEMD-ARIMA 0.431 0.333 446.648 0.826 0.308 0.244 105.330 0.860 
SPI 6 ARIMA 0.708 0.509 179.427 0.554 0.503 0.546 149.074 0.578 
EEMD-ARIMA 0.714 0.505 193.719 0.672 0.211 0.170 119.262 0.926 
SPI 9 ARIMA 0.495 0.339 395.772 0.775 0.330 0.239 179.918 0.775 
EEMD-ARIMA 0.205 0.144 128.230 0.962 0.108 0.080 57.977 0.976 
SPI 12 ARIMA 0.306 0.182 90.516 0.914 0.214 0.128 58.442 0.899 
EEMD-ARIMA 0.118 0.084 85.546 0.987 0.060 0.047 15.188 0.992 
NameModelTraining
Testing
RMSEMAEMAPER2RMSEMAEMAPER2
SPI 3 ARIMA 0.811 0.632 271.170 0.380 0.658 0.510 387.054 0.353 
EEMD-ARIMA 0.431 0.333 446.648 0.826 0.308 0.244 105.330 0.860 
SPI 6 ARIMA 0.708 0.509 179.427 0.554 0.503 0.546 149.074 0.578 
EEMD-ARIMA 0.714 0.505 193.719 0.672 0.211 0.170 119.262 0.926 
SPI 9 ARIMA 0.495 0.339 395.772 0.775 0.330 0.239 179.918 0.775 
EEMD-ARIMA 0.205 0.144 128.230 0.962 0.108 0.080 57.977 0.976 
SPI 12 ARIMA 0.306 0.182 90.516 0.914 0.214 0.128 58.442 0.899 
EEMD-ARIMA 0.118 0.084 85.546 0.987 0.060 0.047 15.188 0.992 

Table 5 gives a thorough breakdown of how well the ARIMA and EEMD-ARIMA models performed throughout both the training and validation phases in predicting SPI values at various temporal scales.

In the context of the training phase for SPI 3, it is noted that the EEMD-ARIMA model outperforms the ARIMA model, with superior values for RMSE (0.431 versus 0.811), MAE (0.333 versus 0.632), and R2 (0.826 versus 0.380). However, it is acknowledged that the MAPE for EEMD-ARIMA is higher than that of ARIMA. Despite this anomaly in MAPE, the other error indices consistently demonstrate the significant improvement of the EEMD-ARIMA model over the conventional ARIMA approach. EEMD-ARIMA continues to perform better than ARIMA throughout the SPI 3 validation phase (RMSE: 0.658, MAE: 0.510, MAPE: 387.54, R2: 0.353 vs. RMSE: 0.308, MAE: 0.244, MAPE: 105.330, R2: 0.860). While also attaining a higher R2, EEMD-ARIMA demonstrated a substantially lower RMSE, MAE, and MAPE, indicating superior accuracy and a better fit to the validation data.

Upon a detailed examination for training data in SPI 6, the comparison between the ARIMA and EEMD-ARIMA models reveals distinct patterns in their performance metrics. In the ARIMA model, we observe a lower MAPE of 179.427, indicating a relatively better accuracy in predicting absolute percentage errors. On the other hand, the EEMD-ARIMA model shows a higher MAPE of 193.719, suggesting a marginally reduced accuracy in this specific measure. Interestingly, when considering R2 values, the EEMD-ARIMA model outshines the ARIMA model significantly, boasting a much higher R2 of 0.672 compared to 0.554 in ARIMA. This underscores the nuanced strengths and weaknesses of each model, with the ARIMA model excelling in MAPE, while the EEMD-ARIMA model demonstrates superior explanatory power as reflected by the higher R2 values. However, EEMD-ARIMA performed significantly better than ARIMA in the validation phase for SPI 6 (RMSE: 0.503, MAE: 0.546, MAPE: 149.074, R2: 0.578 versus RMSE: 0.211, MAE: 0.170, MAPE: 119.262, R2: 0.926). In addition to generating a significantly higher R2, EEMD-ARIMA showed noticeably reduced RMSE, MAE, and MAPE values, indicating enhanced accuracy and a better fit to the validation data.

In the training phase for SPI 9, both models showed relatively good performance. ARIMA (RMSE: 0.495, MAE: 0.339, MAPE: 395.722, R2: 0.775) and EEMD-ARIMA (RMSE: 0.205, MAE: 0.144, MAPE: 128.230, R2: 0.962) exhibited differing RMSE and MAE values, but both had reasonable R2 values, indicating good agreement with the training data. In the validation phase for SPI 9, EEMD-ARIMA (RMSE: 0.108, MAE: 0.080, MAPE: 57.977, R2: 0.976) continued to outperform ARIMA (RMSE: 0.330, MAE: 0.239, MAPE: 179.218, R2: 0.775) by a significant margin. EEMD-ARIMA showed substantially lower RMSE, MAE, and MAPE values, indicating superior accuracy and a better fit to the validation data, while also achieving a considerably higher R2.

For SPI 12, both models performed well during the training phase. In spite of having different RMSE and MAE values, both ARIMA (RMSE: 0.306, MAE: 0.182, MAPE: 90.516, R2: 0.914) and EEMD-ARIMA (RMSE: 0.118, MAE: 0.084, MAPE: 85.546, R2: 0.978) produced high R2 values, indicating an excellent fit to the training data. In the validation phase for SPI 12, EEMD-ARIMA (RMSE: 0.214, MAE: 0.128, MAPE: 58.442, R2: 0.899) considerably beats ARIMA (RMSE: 0.060, MAE: 0.047, MAPE: 15.188, R2: 0.992). In addition to generating a significantly better R2, EEMD-ARIMA demonstrated noticeably reduced RMSE, MAE, and MAPE values, indicating superior accuracy and a superb fit to the validation data.

The result of this analysis reveals that the EEMD-ARIMA model consistently outperformed the ARIMA model across all SPI timescales and in both training and validation phases. EEMD-ARIMA also consistently achieved higher R2 values and showed significantly lower RMSE, MAE, and MAPE values, indicating higher accuracy and better goodness of fit to the data.

Results of SPI predictions using the ARIMA and EEMD-ARIMA models are shown in Figures 811. These figures show that the EEMD-ARIMA model can simulate SPI more accurately than the ARIMA model. This conclusion asserts the suitability of the EEMD model for decomposing monthly precipitation time series, underscores the viability of the ‘decomposition and ensemble’ concept, and highlights the potential of the proposed EEMD-ARIMA model to overcome the limitations of individual models by generating a synergistic effect in forecasting. These claims find support in the presented data, where the EEMD-ARIMA model consistently outperforms the individual ARIMA model across various error indices, demonstrating its enhanced accuracy in predicting drought events. The effectiveness of the decomposition and ensemble approach is evident in the improved forecast precision, as illustrated by the reduced RMSE, MAE, and increased R2 values. Thus, the robustness of the EEMD-ARIMA model is not only a theoretical proposition but is substantiated by the empirical evidence provided through our comprehensive analysis of the monthly precipitation dataset. Thus, using monthly precipitation and corresponding SPI data that have been deconstructed using the EEMD technique as model input data can increase prediction accuracy.
Figure 8

ARIMA, EEMD-ARMA forecasted and observed SPI 3 during validation period.

Figure 8

ARIMA, EEMD-ARMA forecasted and observed SPI 3 during validation period.

Close modal
Figure 9

ARIMA, EEMD-ARMA forecasted and observed SPI 6 during validation period.

Figure 9

ARIMA, EEMD-ARMA forecasted and observed SPI 6 during validation period.

Close modal
Figure 10

ARIMA, EEMD-ARMA forecasted and observed SPI 9 during validation period.

Figure 10

ARIMA, EEMD-ARMA forecasted and observed SPI 9 during validation period.

Close modal
Figure 11

ARIMA, EEMD-ARMA forecasted and observed SPI 12 during validation period.

Figure 11

ARIMA, EEMD-ARMA forecasted and observed SPI 12 during validation period.

Close modal

The efficiency of the EEMD preprocessing method is demonstrated by the EEMD-ARIMA model's persistent outperformance over the ARIMA model in SPI-based drought prediction, as shown by much lower RMSE, MAE, and MAPE values together with higher R2 values. Even though there are sporadic higher MAPE occurrences, the general improvement highlights how well EEMD and ARIMA work together. By utilizing EEMD's advantages in managing nonstationary and nonlinear components, the hybrid model allows ARIMA to focus on linear trends and autocorrelations. This cooperative method produces a greater accuracy synergy, which is especially noticeable during the validation stages. The theoretical claim is supported by the empirical data, which highlights the EEMD-ARIMA model's resilience in improving drought forecast accuracy by overcoming the shortcomings of individual models. To sum up, this clearly illustrates the ‘decomposition and ensemble’ notion in SPI-based drought forecasting, highlighting the importance of EEMD incorporation as a preprocessing step to improve overall prediction accuracy.

Due to the requirement for parameter tuning, which includes optimizing EEMD ensemble trials, ARIMA orders, and decomposition levels, EEMD-ARIMA has some computational challenges that stem from its heavy computational load. Addressing missing data and nonstationarities adds to this computational load, too. Because of the randomness introduced by the ensemble component, consistent findings need several runs. There are difficulties in interpreting a potentially large number of decomposed components, and the long training times of ARIMA models for individual components could be a limitation.

In areas with comparable climates, the EEMD-ARIMA approach has the potential to increase the precision of drought forecasts. Its versatility in capturing nonstationary and nonlinear precipitation patterns – shown in Afghanistan's Herat province – indicates that it might find use elsewhere. However, because different climate dynamics exist in different places, its performance may vary. It is important to weigh robustness and constraints, highlighting the requirement for region-specific calibration to guarantee accurate forecasts. Variations in precipitation patterns and environmental factors may have an impact on the model's efficacy. Furthermore, different geographical contexts may demand different resource requirements and computing efficiency, thus deploying the model there will require careful attention.

Modern water resources management has the crucial but frequently challenging issue of increasing the accuracy of hydrological forecasts, particularly long-term drought forecasting. The ARIMA model has emerged as one of the most widely used approaches in forecasting theory and practice as a way of Box–Jenkins time series analysis. In this work, an effort was made to evaluate the effectiveness of the ARIMA method based on EEMD for forecasting drought time series to increase its forecasting accuracy. As a benchmark for comparison, an ARIMA model built on the original SPI time series is also used. The models described in this article were developed using monthly precipitation data from Herat, Afghanistan. To assess the performances of the models produced, four common statistical performance evaluation metrics (RMSE, MAE, MAPE, and R2) are used.

The findings of this study show that the suggested hybrid EEMD-ARIMA model can greatly improve ARIMA time series techniques for drought prediction. The suggested approach has a number of benefits. First, although EEMD's fundamental idea is straightforward, it can offer profound insight into the properties of time series representing SPI. Second, only data from the relevant SPI time series are needed for ARIMA forecasting. Third, conducting ARIMA modeling benefits from knowing the zero mean of IMF components. Last but not least, the suggested models do not necessitate laborious decision-making on the explicit form of models for each unique case. As a result, using a hybrid forecasting model that includes EEMD presents findings that are more accurate and reliable. It may also be useful in studies that focus on hydrological time series forecasting for a variety of issues connected to successful reservoir management. In comparison to earlier studies, our study marks a substantial leap in drought forecasting, extending the application of the EEMD-ARIMA model to Herat's particular environment and giving a blueprint for similar places confronting analogous issues. Our work illustrates the possibility of novel methodology that unite conventional hydrological modeling with modern data-driven techniques, enhancing the conversation on climate adaption options by using the EEMD and the SPI time series.

The first author would like to express gratitude to Afghanistan's Ministry of Higher Education (MoHE) for the scholarship and Kabul Education University (KEU) for the study leave.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.

Afghanistan: Drought - Emergency Plan of Action (EPoA) DREF Operation n° MDRAF007. 2021 International Federation of Red Cross and Red Crescent Societies.
https://reliefweb.int/report/afghanistan/afghanistan-drought-emergency-plan-action-epoa-dref-operation-n-mdraf007.
Afghanistan: Drought response. 2018 Food and Agriculture Organization of the United Nations.
https://reliefweb.int/report/afghanistan/afghanistan-drought-response-november-2018.
Akaike
H.
1969
Fitting autoregressive models for prediction
.
Ann. Inst. Stat. Math.
21
(
1
),
243
247
.
doi: 10.1007/BF02532251
.
Alam
N. M.
,
Mishra
P.
&
Adhikary
P. P.
2014
Stochastic model for drought forecasting for Bundelkhand region in Central India mapping and characterization of Jhola land (modified stream bed) areas in Koraput district (Eastern Ghats High Land Region of Odisha) view project modelling groundwater nitrate pollution at IARI farm view project. Available at https://www.researchgate.net/publication/261136160.
Alawsi
M. A.
,
Zubaidi
S. L.
,
Al-Bdairi
N. S. S.
,
Al-Ansari
N.
&
Hashim
K.
2022
Drought forecasting: A review and assessment of the hybrid techniques and data pre-processing
.
Hydrology
9
(
7
),
115
.
doi: 10.3390/hydrology9070115
.
Anshuka
A.
,
van Ogtrop
F. F.
&
Willem Vervoort
R.
2019
Drought forecasting through statistical models using standardised precipitation index: A systematic review and meta-regression analysis
.
Nat. Hazards
97
(
2
).
doi: 10.1007/s11069-019-03665-6
.
Badrzadeh
H.
,
Sarukkalige
R.
&
Jayawardena
A. W.
2016
Improving ANN-based short-term and long-term seasonal river flow forecasting with signal processing techniques
.
River Res. Appl.
32
(
3
),
245
256
.
doi: 10.1002/rra.2865
.
Başakın
E. E.
,
Ekmekcioğlu
Ö.
&
Özger
M.
2021
Drought prediction using hybrid soft-computing methods for semi-arid region
.
Model. Earth Syst. Environ.
7
(
4
),
2363
2371
.
doi: 10.1007/s40808-020-01010-6
.
Belayneh
A.
&
Adamowski
J.
2012
Standard precipitation index drought forecasting using neural networks, wavelet neural networks, and support vector regression
.
Appl. Comput. Intell. Soft Comput.
2012
,
1
13
.
doi: 10.1155/2012/794061
.
Box
G. E.
&
Jenkins
G.
1976
Time Series Analysis: Forecasting and Control
.
Holden-Day
,
San Francisco, CA, USA
.
Box
G. E. P.
,
Jenkins
G. M.
&
Reinsel
G. C.
1994
Time Series Analysis: Forecasting and Control
, 3rd edn.
Prentice Hall
,
New Jersey
, USA.
Buyukyildiz
M.
,
Tezel
G.
&
Yilmaz
V.
2014
Estimation of the change in lake water level by artificial intelligence methods
.
Water Resour. Manag.
28
(
13
),
4747
4763
.
doi: 10.1007/s11269-014-0773-1
.
Cacciamani
C.
,
Morgillo
A.
,
Marchesi
S.
&
Pavan
V.
2007
Monitoring and forecasting drought on a regional scale: Emilia-Romagna region
. In: Rossi, G., Vega, T. & Bonaccorso, B. (eds).
Methods and Tools for Drought Analysis and Management
.
Springer Netherlands
,
Dordrecht
, Netherlands pp.
29
48
.
Cancelliere
A.
,
Di Mauro
G.
,
Bonaccorso
B.
&
Rossi
G.
2007
Drought forecasting using the standardized precipitation index
.
Water Resour. Manag.
21
(
5
).
doi: 10.1007/s11269-006-9062-y
.
Chan
W.-S.
1999
A comparison of some of pattern identification methods for order determination of mixed ARMA models
.
Stat. Probab. Lett.
42
(
1
),
69
79
.
doi: 10.1016/S0167-7152(98)00195-3
.
Das
P.
,
Naganna
S. R.
,
Deka
P. C.
&
Pushparaj
J.
2020
Hybrid wavelet packet machine learning approaches for drought modeling
.
Environ. Earth Sci.
79
(
10
),
221
.
doi: 10.1007/s12665-020-08971-y
.
Di
C.
,
Yang
X.
&
Wang
X.
2014
A four-stage hybrid model for hydrological time series forecasting
.
PLoS One
9
(
8
),
e104663
.
doi: 10.1371/journal.pone.0104663
.
Djerbouai
S.
&
Souag-Gamane
D.
2016
Drought forecasting using neural networks, wavelet neural networks, and stochastic models: case of the Algerois basin in North Algeria
.
Water Resour. Manag.
30
(
7
),
2445
2464
.
doi: 10.1007/s11269-016-1298-6
.
Elliott
G.
,
Rothenberg
T. J.
&
Stock
J. H.
1996
Efficient tests for an autoregressive unit root
.
Econometrica
64
(
4
),
813
.
doi: 10.2307/2171846
.
Gumus
V.
&
Algin
H. M.
2017
Meteorological and hydrological drought analysis of the Seyhan − Ceyhan River Basins, Turkey
.
Meteorol. Appl.
24
(
1
),
62
73
.
doi: 10.1002/met.1605
.
Guttman
N. B.
1998
Comparing the palmer drought index and the standardized precipitation index
.
J. Am. Water Resour. Assoc.
34
(
1
).
doi: 10.1111/j.1752-1688.1998.tb05964.x
.
Guttman
N. B.
1999
Accepting the standardized precipitation index: A calculation algorithm
.
JAWRA J. Am. Water Resour. Assoc.
35
(
2
),
311
322
.
doi: 10.1111/j.1752-1688.1999.tb03592.x
.
Hayes, M., Svoboda, M., Le Comte, D., Redmond, K. T. & Pasteris, P. 2005 Drought monitoring: New tools for the 21st century. Drought and Water Crises: Science, technology, and management issues 53, 69.
Hayes
M.
,
Svoboda
M.
,
Wall
N.
&
Widhalm
M.
2011
The Lincoln declaration on drought indices: Universal meteorological drought index recommended
.
Bull. Am. Meteorol. Soc.
92
(
4
).
doi: 10.1175/2010BAMS3103.1
.
Hou
H.
&
Andrews
H.
1978
Cubic splines for image interpolation and digital filtering
.
IEEE Trans. Acoust.
26
(
6
),
508
517
.
doi: 10.1109/TASSP.1978.1163154
.
Huang
N. E.
&
Wu
Z.
2008
A review on Hilbert-Huang transform: Method and its applications to geophysical studies
.
Rev. Geophys.
46
(
2
),
RG2006
.
doi: 10.1029/2007RG000228
.
Huang
N. E.
,
Shen
Z.
,
Long
S. R.
,
Wu
M. C.
,
Shih
H. H.
,
Zheng
Q.
,
Yen
N. C.
,
Tung
C. C.
&
Liu
H. H.
1998
The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis
.
Proc. R. Soc. London. Ser. A Math. Phys. Eng. Sci.
454
(
1971
),
903
995
.
doi: 10.1098/rspa.1998.0193
.
Khan
M. M. H.
,
Muhammad
N. S.
&
El-Shafie
A.
2020
Wavelet based hybrid ANN-ARIMA models for meteorological drought forecasting
.
J. Hydrol.
590
,
125380
.
doi: 10.1016/j.jhydrol.2020.125380
.
Kisi
O.
,
Latifoğlu
L.
&
Latifoğlu
F.
2014
Investigation of empirical mode decomposition in forecasting of hydrological time series
.
Water Resour. Manag.
28
(
12
),
4045
4057
.
doi: 10.1007/s11269-014-0726-8
.
Koycegiz
C.
&
Buyukyildiz
M.
2022
Investigation of precipitation and extreme indices spatiotemporal variability in Seyhan Basin, Turkey
.
Water Supply
22
(
12
),
8603
8624
.
doi: 10.2166/ws.2022.391
.
Kumar
R.
, Musuuza, J. L., Van Loon, A. F., Teuling, A. J., Barthel, R., Ten Broek, J., Mai, J., Samaniego, L. & Attinger, S.
2016
Multiscale evaluation of the standardized precipitation index as a groundwater drought indicator
.
Hydrol. Earth Syst. Sci.
20
(
3
),
1117
1131
.
doi: 10.5194/hess-20-1117-2016
.
Lee
T.
&
Ouarda
T. B. M. J.
2010
Long-term prediction of precipitation and hydrologic extremes with nonstationary oscillation processes
.
J. Geophys. Res.
115
(
D13
),
D13107
.
doi: 10.1029/2009JD012801
.
Libanda
B.
&
Nkolola
N. B.
2022
An ensemble empirical mode decomposition of consecutive dry days in the Zambezi Riparian Region: Implications for water management
.
Phys. Chem. Earth, Parts A/B/C
126
,
103147
.
doi: 10.1016/j.pce.2022.103147
.
Liu
Y.
,
Brown
J.
,
Demargne
J.
&
Seo
D.-J.
2011
A wavelet-based approach to assessing timing errors in hydrologic predictions
.
J. Hydrol.
397
(
3–4
),
210
224
.
doi: 10.1016/j.jhydrol.2010.11.040
.
Liu
D.
,
You
J.
,
Xie
Q.
,
Huang
Y.
&
Tong
H.
2018
Spatial and temporal characteristics of drought and flood in Quanzhou based on standardized precipitation index (SPI) in recent 55 years
.
J. Geosci. Environ. Prot.
06
(
08
),
25
37
.
doi: 10.4236/gep.2018.68003
.
Liu
X.
,
Zhang
Y.
&
Zhang
Q.
2022
Comparison of EEMD-ARIMA, EEMD-BP and EEMD-SVM algorithms for predicting the hourly urban water consumption
.
J. Hydroinformatics
24
(
3
),
535
558
.
doi: 10.2166/hydro.2022.146
.
Mayar
M. A.
2021
Droughts on the Horizon: Can Afghanistan manage this risk? Afghanistan-analysts.org
.
McKee
T. B.
,
Doesken
N.
&
Kleist
J.
1993
The relationship of drought frequency and duration to time scales. Proceedings of the 8th Conference on Applied Climatology 17 (22), 179–183
.
Mishra
A. K.
&
Desai
V. R.
2005
Drought forecasting using stochastic models
.
Stoch. Environ. Res. Risk Assess.
19
(
5
),
326
339
.
doi: 10.1007/s00477-005-0238-4
.
Moriasi
D. N.
,
Arnold
J. G.
,
Van Liew
M. W.
,
Bingner
R. L.
,
Harmel
R. D.
&
Veith
T. L.
2007
Model evaluation guidelines for systematic quantification of accuracy in watershed simulations
.
Trans. ASABE
50
(
3
),
885
900
.
doi: 10.13031/2013.23153
.
Nourani
V.
,
Alami
M. T.
&
Aminfar
M. H.
2009
A combined neural-wavelet model for prediction of Ligvanchai watershed precipitation
.
Eng. Appl. Artif. Intell.
22
(
3
),
466
472
.
doi: 10.1016/j.engappai.2008.09.003
.
Nourani
V.
,
Komasi
M.
&
Alami
M. T.
2012
Hybrid wavelet–genetic programming approach to optimize ANN modeling of rainfall–runoff process
.
J. Hydrol. Eng.
17
(
6
),
724
741
.
doi: 10.1061/(ASCE)HE.1943-5584.0000506
.
Nourani
V.
,
Hosseini Baghanam
A.
,
Adamowski
J.
&
Kisi
O.
2014
Applications of hybrid wavelet–artificial intelligence models in hydrology: A review
.
J. Hydrol.
514
,
358
377
.
doi: 10.1016/j.jhydrol.2014.03.057
.
Ntale
H. K.
&
Gan
T. Y.
2003
Drought indices and their application to East Africa
.
Int. J. Climatol.
23
(
11
),
1335
1357
.
doi: 10.1002/joc.931
.
Pham
Q. B.
,
Yang
T.-C.
,
Kuo
C.-M.
,
Tseng
H.-W.
&
Yu
P.-S.
2021
Coupling singular spectrum analysis with least square support vector machine to improve accuracy of SPI drought forecasting
.
Water Resour. Manag.
35
(
3
),
847
868
.
doi: 10.1007/s11269-020-02746-7
.
Prajapati, V. K., Khanna, M., Singh, M., Kaur, R., Sahoo, R. N., Singh, D. K. & Datta, S. P. 2021 Drought monitoring using multi-time-scale Standardized Precipitation Index. The Indian Journal of Agricultural Sciences 91 (10), 1438–1442. https://doi.org/10.56093/ijas.v91i10.117423.
Rezaiy
R.
&
Shabri
A.
2023
Drought forecasting using W-ARIMA model with standardized precipitation index
.
J. Water Clim. Chang.
14
(
9
),
3345
3367
.
doi: 10.2166/wcc.2023.431
.
Rossi, G. 2003 Requisites for a drought watch system. Tools for Drought Mitigation in Mediterranean Regions 44, 147–157
.
Ruleman
C. A.
,
Crone
A. J.
,
Machette
M. N.
,
Haller
K. M.
&
Rukstales
K. S.
2007
Map and database of probable and possible quaternary faults in Afghanistan
. US Geological Survey open-file report, 1103 (1). USGS Open-File Report 2007-1103: Map and Database of Probable and Possible Quaternary Faults in Afghanistan.
Schwarz
G.
1978
Estimating the dimension of a model
.
Ann. Stat.
6
(
2
),
461
464
.
Shabri
A.
2014
A hybrid wavelet analysis and adaptive neuro-fuzzy inference system for drought forecasting
.
Appl. Math. Sci.
8
,
6909
6918
.
doi: 10.12988/ams.2014.48263
.
Sönmez
F. K.
,
Kömüscü
A. Ü.
,
Erkan
A.
&
Turgu
E.
2005
An analysis of spatial and temporal dimension of drought vulnerability in Turkey using the standardized precipitation index
.
Nat. Hazards
35
(
2
),
243
264
.
doi: 10.1007/s11069-004-5704-7
.
Timuri
N. a-D.
2023
Investigating drought in herat province and ways to combat It
.
Jami Sci. Res. Q. J.
8
(
2
),
55
66
.
doi: 10.61438/jsrqj.v8i2.27
.
Tongal
H.
2013
Nonlinear forecasting of stream flows using a chaotic approach and artificial neural networks
.
Earth Sci. Res. J.
17
(
2
),
119
126
.
Wang
W.
,
Chau
K.
,
Xu
D.
&
Chen
X.-Y.
2015
Improving forecasting accuracy of annual runoff time series using ARIMA based on EEMD decomposition
.
Water Resour. Manag.
29
(
8
),
2655
2675
.
doi: 10.1007/s11269-015-0962-6
.
Wang
Z.-Y.
,
Qiu
J.
&
Li
F.-F.
2018
Hybrid models combining EMD/EEMD and ARIMA for long-term streamflow forecasting
.
Water
10
(
7
),
853
.
doi: 10.3390/w10070853
.
Wu
Z.
&
Huang
N. E.
2009
Ensemble empirical mode decomposition: a noise-assisted data analysis method
.
Adv. Adapt. Data Anal.
01
(
01
),
1
41
.
doi: 10.1142/S1793536909000047
.
Wu
X.
,
Zhou, J., Yu, H., Liu, D., Xie, K., Chen, Y., Hu, J., Sun, H. & Xing, F.
2021
The development of a hybrid wavelet-ARIMA-LSTM model for precipitation amounts and drought analysis
.
Atmosphere (Basel)
12
(
1
),
74
.
doi: 10.3390/atmos12010074
.
Xu
D.
,
Ding
Y.
,
Liu
H.
,
Zhang
Q.
&
Zhang
D.
2022
Applicability of a CEEMD–ARIMA combined model for drought forecasting: A case study in the Ningxia Hui Autonomous Region
.
Atmosphere (Basel)
13
(
7
),
1109
.
doi: 10.3390/atmos13071109
.
Yihdego
Y.
,
Vaheddoost
B.
&
Al-Weshah
R. A.
2019
Drought indices and indicators revisited
.
Arab. J. Geosci.
12
(
3
).
doi: 10.1007/s12517-019-4237-z
.
Zargar
A.
,
Sadiq
R.
,
Naser
B.
&
Khan
F. I.
2011
A review of drought indices
.
Environ. Rev.
19
.
doi: 10.1139/a11-013
.
Zhang
P. G.
2003
Time series forecasting using a hybrid ARIMA and neural network model
.
Neurocomputing
50
,
159
175
.
doi: 10.1016/S0925-2312(01)00702-0
.
Zhang
J.
,
Yan
R.
,
Gao
R. X.
&
Feng
Z.
2010
Performance enhancement of ensemble empirical mode decomposition
.
Mech. Syst. Signal Process.
24
(
7
),
2104
2123
.
doi: 10.1016/j.ymssp.2010.03.003
.
Zhang
H.
,
Singh
V. P.
,
Wang
B.
&
Yu
Y.
2016
CEREF: A hybrid data-driven model for forecasting annual streamflow from a socio-hydrological system
.
J. Hydrol.
540
,
246
256
.
doi: 10.1016/j.jhydrol.2016.06.029
.
Zhao
X.
&
Chen
X.
2015
Auto regressive and ensemble empirical mode decomposition hybrid model for annual runoff forecasting
.
Water Resour. Manag.
29
(
8
),
2913
2926
.
doi: 10.1007/s11269-015-0977-z
.
This is an Open Access article distributed under the terms of the Creative Commons Attribution Licence (CC BY 4.0), which permits copying, adaptation and redistribution, provided the original work is properly cited (http://creativecommons.org/licenses/by/4.0/).