Forecasting of Standardized Precipitation Index using Hybrid Models: A Case Study of Cape Town, South Africa

1. Introduction

Drought is a gradual and pernicious natural catastrophic event with global socioeconomic and environmental consequences. This is a highly perilous climate-related catastrophe that has a substantial effect on both the environment and human existence [1]. Droughts have a profound effect on municipal water supplies, which is one of their most important effects. These supplies refer to the purified water distributed to households and businesses, which is crucial for human consumption, cleanliness, and numerous critical services on a global scale [2]. Municipal water sources are given priority during drought and water constraint due to their significant sociological and economic relevance. This prioritization entails reallocating available water resources, frequently diverting them away from other sectors, especially agriculture, to ensure the uninterrupted supply of water to urban areas and their residents [3]. During extremely severe drought circumstances, the water sources that urban water systems depend on can be exhausted, leading to a scarcity of municipal water for consumers.

In 2018, Cape Town and another part of the country faced acute water shortages as a result of drought. These occurrences and growing concerns about water scarcity in urban areas are the result of several factors, including rising demand for water due to population growth, unpredictable weather patterns, the effects of climate change, and, in some cases, insufficient planning and management of water infrastructure [4,5]. In addition, the exploitation of surface water and groundwater resources, which have historically supplied water for human sustenance and agriculture, has reached unsustainable levels in several countries, resulting in negative environmental consequences [6]. The excessive use of these natural water resources has led to the contamination of certain water bodies, such as the salinization of over-exploited aquifers or the inability to dilute and assimilate wastewater discharges, rendering the water supplies unsuitable for human consumption [7]. To maintain a secure and stable water supply in the present and the future, it is essential to recognize that the development of resilient water systems is an imperative necessity [8]. In disaster management, the application of time series forecasting methodologies can serve as an early warning system.

The significance of addressing drought-related issues has increased as time series analysis continues to advance. Prediction of drought requires the management of complex, often nonstationary, nonlinear precipitation data, and SPI time series. Rarely do we encounter simple and stable time series in practice. Consequently, it is of the utmost importance to identify an efficient method for addressing the complexities of nonlinear and nonstationary time series and for making accurate forecasts. Nonetheless, due to the complexities of evaluating time series precision, devising suitable forecasting models remains a formidable challenge [9].

In recent years, numerous models, such as the Autoregressive Integrated Moving Average (ARIMA) model, have been utilized for drought forecasting due to its their adaptability and greater information on time-related changes [10]. However, they have difficulty accurately predicting nonlinear data [11]. As machine learning has evolved, the Long Short-Term Memory (LSTM) network has emerged as a solution to resolve long-term dependencies, especially in sequences with extended intervals and delays [12]. The LSTM model has shown promising results when applied to the prediction of drought time series and related domains [13,14,15]. Nonetheless, the inherent complexity of time series frequently results in suboptimal local predictions, thereby diminishing the overall performance of forecasting [11].

To enhance time-series prediction accuracy in the face of these complexities, scholars have introduced signal decomposition techniques [16,17,18,19]. Signal decomposition efficiently removes certain characteristics from sequences, making them more stable and easier to handle. These methods streamline the initial simplify time series and enhance their predictability. Some of the methods include empirical mode decomposition (EMD), ensemble empirical mode decomposition (EEMD), complementary ensemble empirical mode decomposition (CEEMD), and complementary ensemble empirical mode decomposition with adaptive noise (CEEMDAN). CEEMD addresses both the problem of blending modes in EMD and the issue of residual white noise in EEMD. Furthermore, CEEMDAN is an enhanced iteration of CEEMD that provides enhanced adaptability in handling intricate signals with varying attributes and levels of noise. CEEMDAN is specifically intended to adaptively assess and remove noise from the signal throughout the decomposition process. This is especially advantageous when the input signal is contaminated by noise, as it enhances the integrity of the recovered components. Xu et al. [20] utilised the CEEMD-ARIMA combined model to forecast drought in the Xinjiang Uygur Autonomous Region of China. The model outperformed the standalone ARIMA in terms of prediction accuracy across many time scales. Ding et al. [11] utilised a hybrid model combining CEEMD-LSTM to predict drought in the Xinjiang Uygur Autonomous Region of China using the standardised precipitation index. The hybrid model demonstrated superior performance compared to LSTM.

Given the non-linear and non-stationary characteristics of precipitation data, we investigate more efficient approaches for predicting time-series. Within this framework, we suggest a hybrid model called CEEMDAN-ARIMA-LSTM. This model combines the integrated empirical mode decomposition with adaptive noise (CEEMDAN), the differential integrated moving autoregressive average model (ARIMA), and the Long Short-Term Memory network (LSTM). The objective is to enhance the precision of drought prediction by capitalising on the strengths of these methods in managing intricate time series. By conducting a comparison analysis between the LSTM and ARIMA approaches, we have determined that our combined model is capable of successfully analysing and controlling the specific features of time series. This allows us to utilise both explicit and hidden data information, resulting in a substantial improvement in prediction accuracy.

2. Materials and Methods

2.1. Study Area

The data used in the study is from Cape Town as shown in Figure 1, also known as "Mother City," that which is located in the Western Cape province of South Africa, near the southernmost point of the African continent. The exact geographical coordinates are approximately 33.9249° S latitude and 18.4241° E longitude. This coastal city is located on the southwestern coast of South Africa, between the Atlantic Ocean to the west and the Hottentots-Holland Mountains to the east, creating a breathtaking natural environment that includes picturesque coastlines and Table Mountain. The unique geographical location of Cape Town has played a significant role in its history and development, contributing to the city's breathtaking natural attractiveness.

2.2. Mann-Kendall Test Statistics

Examining potential changes in the characteristics of precipitation events resulting from climate change is of utmost importance as it has significant implications for soil moisture retention, hazard assessment, and flood management. The assessment of monotonous trends in a time series of geophysical data proves consistently advantageous. The Mann-Kendall (MK) test and modified Mann-Kendall (mMK), are widely used nonparametric statistical methods, serving as an effective tool for evaluating trends in time series data. Mann first introduced this test in 1945 [21], and Kendall later refined it in 1975 [22]. The computation of the test statistic $S$ follows the formula below.

$$S=\sum _{j=1}^{n-1}\sum _{i=j+1}^{n}Sign\left(x_i-x_j\right)$$

(1)

$$Sign\left(x\right)=f\left(x\right)=\left\{\begin{array}{c}1 for x>1\\ 0 for x=1\\ -1, for x<1\end{array}\right.$$

(2)

Let n represent the length of the data. The data values at times $i$ and $j$ are denoted as $x_i$ and $x_j$, respectively. According to Seenu and Jayakumar [23], in the Mann-Kendall test, a positive test statistic, S, suggests an increasing trend, whereas a negative test statistic indicates a decreasing trend. The variance for S is calculated using a specific mathematical formula.

$$Var\left(S\right)={\displaystyle \frac{n\left(n-1\right)\left(2n+5\right)-\sum _{i=1}^p{t}_i\left(t_i-1\right)\left(2t_i+5\right)}{18}}$$

(3)

The variable $t_i$ represents the number of data points in the $j^{th}$ tied group, while $p$ represents the number of the tied group in the time series. The statement indicates that the data in a time series can be compared to each other [24]. The summation procedure in the equation (3) is specifically used for linked groups in the time series. Its purpose is to minimize the impact of individual values within tied groups on the ranking statistics [25]. Assuming random and independent time series, the S statistic follows a normal distribution and can be calculated as the standardized Z statistic as follows

$$Z=\left\{\begin{array}{c}{\displaystyle \frac{S-1}{Var\left(S\right)}} , S>0\\ 0 , S=0 \\ {\displaystyle \frac{S+1}{Var\left(S\right)}} , S<0\end{array}\right.$$

(4)

The Z variable follows a standard normal distribution with a variance of one and a mean of zero. The value of the S statistic is associated with the Kendall term $\tau$, which is defined as

$$\tau ={\displaystyle \frac{S}{D}}$$

(5)

where

$$D={\left[{\displaystyle \frac{1}{2}}n\left(n-1\right)-{\displaystyle \frac{1}{2}}\sum _{j=1}^p{t}_j\left(t_j-1\right)\right]}^{{\displaystyle \frac{1}{2}}}{\left[{\displaystyle \frac{1}{2}}n\left(n-1\right)\right]}^{{\displaystyle \frac{1}{2}}}$$

(6)

The test statistic Z is utilized to quantify the significance of trends. In fact, the null hypothesis of the MK test assumes the absence of any trend and is tested against the alternative hypothesis which assumes the presence of a trend at a specific level of significance [26,27]. Another significant outcome of the Mann-Kendall statistics is the Kendall term, which measures the correlation and indicates the strength of the relationship between any two independent variables. For serially correlated data, having a significant lag-1 autocorrelation coefficient, modified Mann-Kendall test using the Hamed and Ramachandra Rao [28] variance correlation approach was employed [29,30]. The positive and negative values of the standardize MK and MMK test statistics ($Z_{MK}$ and $Z_{MMK})$ indicate increasing and decreasing trends, respectively.

The Mann-Kendall trend method can be further extended to a sequential version of the Mann-Kendall test statistic known as the Sequential Mann-Kendall (SQ-MK) test [31] with an aim to detect trends in time series [32]. The test involves comparing the original time series with its reversed counterpart to identify potential trends. The visual representation of the direct series $u( t_i)$ and the backward series $u^{\prime }\left(t_i\right)$ can offer valuable insights into the presence or absence of trends. The Sequential Mann-Kendall test is computed by assigning ranked values, $y_i$, to the initial values in the analysis $\left(x_1,x_2,x_3,\cdots ,x_n\right).$ The magnitudes of $y_i$ $( i=1, 2, 3,...,n )$ are compared with $y_j$ $( j=1, 2, 3,...,i-1 )$. Each comparison results in a count, denoted as $n_i$, for the cases where $y_i>y_j$. Consequently, a statistic $t_i$ can be precisely defined as:

$$t_i=\sum _{j=1}^i{n}_i$$

(7)

The distribution of test statistic $t_i$ has a mean as

$$E\left(t_i\right)={\displaystyle \frac{i\left(i-1\right)}{4}}$$

(8)

and variance as

$$Var\left(t_i\right)={\displaystyle \frac{i\left(i-1\right)\left(2i+5\right)}{72}}$$

(9)

Finaly, the sequential values of a reduced or standardized variable, called statistic $u( t_i)$ is calculated for each of the test statistic variable $t_i$ as follows:

$$u\left( t_i\right)={\displaystyle \frac{t_i-E\left(t_i\right)}{\sqrt{Var\left(t_i\right)}}}$$

(10)

The equation (10) provides the forward sequential statistic, $u\left( t_i\right)$, also known as the progressive statistic, which is estimated using the original time series. The backward sequential statistic, $u^{\prime }\left(t_i\right)$, also known as the retrograde statistic, is calculated in the same way but starting from the end of the series. To estimate $u^{\prime }\left(t_i\right)$, the time series is rearranged so that the last value of the original time series becomes the first $\left(x_1,x_2,x_3,\cdots ,x_n\right)$. The sequential version of the Mann-Kendall test statistic allows for the detection of an approximate beginning of a developing trend. By plotting the $u( t_i)$and $u^{\prime }\left(t_i\right)$ curves, the intersection of these curves can indicate an approximate potential trend turning point. If the intersection of $u( t_i)$and $u^{\prime }\left(t_i\right)$ occurs beyond the 5% level of the standardized statistic, a detectable change at that point in the time series can be inferred. Additionally, if at least one value of the reduced variable is greater than a chosen level of significance of the Gaussian distribution, the null hypothesis is rejected. To enhance the detection and understanding of emerging trend, innovative trend analysis is applied in this study.

2.2. Innovative Trend Analysis

The innovative trend analysis (ITA) is a new technique proposed by Sen in 2012 [33]. The method assigns low, medium, and high values to a time series, and compared to other non-parametric tests, the ITA technique has wide applicability, regardless of distribution assumptions, dataset size, and serial correlation [23]. Due to these advantages, the ITA technique has been extensively utilized in detecting trends in meteorological and hydrological variables. In this method, the first half of the time series is plotted on the x-axis against the second half on the y-axis. It is evident that the monotone increasing (decreasing) trend in the given time series lies above (below) the 1:1 line. In this study, this approach is used to quantify the results of the MK test. The ITA method is a convenient way to extract trend information from available rainfall time series data.

2.3. Standardized Precipitation Index Calculation

The Standard Precipitation Index (SPI) was created by McKee et al. [34] as a probabilistic indicator (see Table 1). The calculation of SPI involves fitting the precipitation data to a distribution function and determining the magnitude based on the quantile point of the normal distribution, which is defined by the cumulative likelihood of a specific amount of precipitation. One primary benefit of the SPI is that it relies exclusively on precipitation data as its input. This feature makes it well-suited for areas with restricted data accessibility [35]. To calculate precise distribution characteristics, the SPI necessitates a specific quantity of precipitation samples, typically exceeding 20 to 30 years of data [36]. The Gamma distribution is typically the most suitable choice for representing observed precipitation data in the process of computing the SPI. The equation provided represents the probability density function of the Gamma distribution.

$$g\left(x\right)={\displaystyle \frac{1}{{\beta }^{\alpha }\mathsf{\Gamma }\left(\alpha \right)}}{x}^{\alpha -1}{e}^{-{\displaystyle \raisebox{1ex}{$x$}\!\left/ \!\raisebox{-1ex}{$\beta $}\right.}}, \mathrm{f}\mathrm{o}\mathrm{r} x>0$$

(11)

where $\alpha >0$ is the shape parameter, $\beta >0$ is the scale parameter, and $x>0$ is the amount of precipitation. $\mathsf{\Gamma }\left(\alpha \right)$ is the value of the Gamma function, a well-known mathematical function defined by the integral:

$$\mathsf{\Gamma }\left(\alpha \right)=\underset{0}{\overset{\infty }{\int }}{y}^{\alpha -1}{e}^{-y}dy$$

(12)

and the parameter estimates of $\alpha$ and $\beta$ with $n$ observations are defined as:

$$\widehat{\alpha }={\displaystyle \frac{1}{4A}}\left(1+\sqrt{1+{\displaystyle \frac{4A}{3}}}\right) \mathrm{a}\mathrm{n}\mathrm{d} \widehat{\beta }={\displaystyle \frac{\overline{x}}{n}}$$

(13)

$$A=\ln \left(\overline{x}\right)-{\displaystyle \frac{\sum \ln \left(x\right)}{n}}$$

(14)

Once the coefficients $\alpha$ and $\beta$ have been estimated, the probability density function $g\left(x\right)$ is integrated with respect to x. The result is an equation that represents the total probability, denoted as G(x), of observing a certain amount of rainfall within a specific month and time interval.

$$G\left(x\right)=\underset{0}{\overset{x}{\int }}g\left(x\right)dx$$

(15)

In arid and semi-arid regions, precipitation is often zero, leading to the following cumulative probability equation:

$$H\left(x\right)=q+\left(1-q\right)G\left(x\right)$$

(16)

where $q={\displaystyle \frac{m}{n}}$, while $m$ represents the number of zero values in precipitation and $n$ represents the number of observations. Lloyd Hughes and Saunders [37] introduced the SPI calculation as follows:

$$SPI=\left\{\begin{array}{c}-\left(t-{\displaystyle \frac{c_0+c_{1}t+c_2{t}^2}{1+d_{1}t+d_2{t}^2+d_3{t}^3}}\right) ,\hspace{1em}\hspace{1em}\hspace{1em}0<H(x)\le 0.5\\ +\left(t-{\displaystyle \frac{c_0+c_{1}t+c_2{t}^2}{1+d_{1}t+d_2{t}^2+d_3{t}^3}}\right) ,\hspace{1em}\hspace{1em}\hspace{1em}0<H(x)\le 1\end{array}\right.$$

(17)

where $x$ is precipitation, $H\left(x\right)$ is the cumulative probability of precipitation observed, and $c_0$, $c_1$, $c_2$, $d_1$, $d_2$, $d_3$ are constants with the following values:

$c_0=2.515517,c_1=0.802853,c_2=0.010328,$

$d_1=1.432788,d_2=0.189269,d_3=0.001308.$

In addition, $t$ is obtained from next mentioned Equation (8):

$$t=\left\{\begin{array}{c}\sqrt{\ln \left({\left({\displaystyle \frac{1}{H\left(x\right)}}\right)}^2\right)} ,\hspace{1em}\hspace{1em}\hspace{1em}0<H(x)\le 0.5\\ \sqrt{\ln \left({\left({\displaystyle \frac{1}{1-H\left(x\right)}}\right)}^2\right)} ,\hspace{1em}\hspace{1em}0<H(x)\le 1\end{array}\right.$$

(18)

 2.2.2. The Complete Ensemble Empirical Mode with Adaptive Noise (CEEMDAN) 

The CEEMDAN is a signal processing technique that was created by Torres et al. [38]. This method breaks down the original signal into intrinsic mode functions (IMFs). The CEEMDAN algorithm incorporates a fixed quantity of adaptable white noise into each decomposition iteration to resolve the issue of modal mixing in EMD decomposition. Additionally, it tackles the problem of incompleteness resulting from EEMD decomposition and improves computational efficiency by increasing the average number of iterations. Set $x\left(t\right)$ as the initial signal and add Gaussian white noise ${\epsilon }_0\omega {\left(t\right)}^{\left(i\right)}$ assuming a normal distribution with mean of 0 and variance of 1. The representation of the i-th sequence expression is as follows:

$$x{\left(t\right)}^{\left(i\right)}=x\left(t\right)+{\epsilon }_0\omega {\left(t\right)}^{\left(i\right)}, i=1, 2, \cdots ,N$$

(19)

where $\omega {\left(t\right)}^{\left(i\right)}$ represents a set of Gaussian white noise sequences with mean of 0 and variance of 1, ${\epsilon }_0$ denotes the standard deviation of Gaussian white noise. The process of CEEMDAN can be given as follows:

• The EMD decomposition is performed for $x{\left(t\right)}^{\left(i\right)}$ to obtain N new sequences and calculate the mean worth to be the first model component IMF₁. The first remaining component $r_1$ will be calculated

$${IMF}_1={\displaystyle \frac{1}{N}}\sum _{i=1}^N{E}_1\left[x{\left(t\right)}^{\left(i\right)}\right]$$

(20)

$$r_1=x\left(t\right)-{IMF}_1$$

(21)

where E is the EMD decomposition operator, $r_1$ denotes the residual signal after the first decomposition.

2.

Adding specific noise to the new signal, EMD decomposition is continued to obtain the second IMF₂ of the original signal and corresponding residual $r_2$

$${IMF}_2={\displaystyle \frac{1}{N}}\sum _{i=1}^N{E}_1\left(r_1+{\epsilon }_1{E}_1\left[\omega {\left(t\right)}^{\left(i\right)}\right]\right)$$

(22)

$$r_2=r_1-{IMF}_2$$

(23)

3.

In the following stage, for $i=3, 4,\cdots ,k$ , the k-th mode component and the corresponding residual signal can be computed in the following equation.

$${IMF}_k={\displaystyle \frac{1}{N}}\sum _{i=1}^N{E}_1\left(r_{k-1}+{\epsilon }_{k-1}{E}_{k-1}\left[\omega {\left(t\right)}^{\left(i\right)}\right]\right)$$

(24)

$$r_k=r_{k-1}-{IMF}_k$$

(25)

4.

Repeat step 3 until the residual satisfies the stoppage criterion.

5.

Finally, the decomposition consequence of the original signal can be described as

$$x\left(t\right)=\sum _{k=1}^K{IMF}_k+r_k$$

(26)

2.4. Autoregressive Integrated Moving Average

ARIMA is a model of mathematical statistics that was built in the early 1970s by Box and Jenkins [39]. It is a well-known method for predicting time series that has attracted a great deal of interest in numerous fields and scientific research. The ARIMA model, which stands for Autoregressive Integrated Moving Average, is an advanced technique used for time-series analysis. It is an extension of the Autoregressive (AR) model and the Moving Average (MA) model. The ARMA model, which is a key component of the ARIMA model, is expressed as follows:

$$X_t=Z_t+\sum _{i=1}^q{\theta }_i{Z}_{t-i}+\sum _{i=1}^p{\varphi }_i{X}_{t-i}$$

(27)

where, $X_t$ is the ARMA process, denoted as $ARMA(p,q);$ ${\varphi }_i$ is the AR parameter; ${\theta }_i$ is the MA parameter; $Z_t ~N(0, {\sigma }_z^2)$ is a Gaussian white nose. The ARMA model can be extended to the ARIMA model by smoothing a nonstationary time series with d-order difference

$${△}^d{X}_t={△}^{d-1}{X}_t-{△}^{d-1}{X}_{t-1}$$

(28)

where $d$ is the order of difference; ${△}^d$ is the $d$-order difference operator. The predictive efficacy of the ARIMA model is contingent upon its parameters $({\varphi }_i, i=\mathrm{1,2},\cdots ,p)$ and (${\theta }_i, i=\mathrm{1,2},\cdots ,q)$, as well as its orders (p, d, and q). The model effectively represents the linear relationship in the series when the parameters and ordering are appropriately specified [40]. The configuration of the ARIMA model is as follows:

• Stationarity test: The augmented Dickey-Fuller (ADF) test is employed to determine whether a time series is stationary or nonstationary. If the time series is nonstationary, it should be subjected to differential operation.
• Order identification: The orders of the ARIMA model are determined by examining the autocorrelation coefficient function (ACF) and partial autocorrelation coefficient function (PACF) of the time series. Another approach Python package is the auto_arima function that conducts a systematic search across potential ARIMA hyperparameters in a stepwise manner.
• Model fitting: The Kalman filter and the maximum likelihood function are used to estimate and fit the model's unknown parameters. Furthermore, the model is examined to see if the residual sequence contains white noise. Otherwise, the orders must be inappropriate and must be re-identified.
• Application: The ARIMA model is utilized to forecast the SPI index, and the data is added to the training sample. Before the next forecast, the model should be updated. Figure 3 depicts the forecasting process using the best ARIMA model.

2.5. Long Short-Term Memory Neural Network

The Long Short-Term Memory (LSTM) was first proposed by Hochreiter and Schmidhuber [41] as one of the most successful recurrent neural networks (RNN) architectures. LSTM, a variant of recurrent neural networks (RNNs), incorporates a threshold mechanism consisting of an input gate, forgetting gate, and output gate. This mechanism evaluates incoming information to determine if it satisfies certain criteria, thus regulating the rate at which information is accumulated. Consequently, LSTM is capable of retaining and updating new information, effectively addressing the issue of long-term dependence. As illustrated in Figure 4, the neural unit of each LSTM is made up of the cell state, namely the long-term state $c_t$ and the short-term state $h_t$, as well as the input gate $i_t$, forgetting gate $f_t$, and output gate $o_t$.

The cell state, which can be viewed as a container for storing information, modifies and outputs the information progressively through process control of the input, forgetting, and output gates. Each cell state within a neural unit experiences the forget gate, the input gate, and the process of transferring information to the output gate. To process the current neural unit, the input gate duplicates the incoming data. The entire input consists of two components. The sigmoid activation function portion specifies which categories of input data are updated, i.e. which input data is disregarded. The hyperbolic tangent function (tanh) was utilised to produce a novel candidate value vector, which was then added to the current cell state. The mathematical procedure is defined by:

$$i_t=\sigma \left(W_i .\left[x_t.h_{t-1}\right]+b_i\right)$$

(29)

$$c_t=\phi \left(W_c . \left[x_t.h_{t-1}\right]+b_c\right)$$

(30)

The forget gate is employed to determine which information in the current state should be discarded, while the LSTM algorithm acquires the ability to instruct the network to retain information.

$$f_t=\sigma \left(W_f .\left[x_t.h_{t-1}\right]+b_f\right)$$

(31)

The output gate is primarily responsible for controlling the output information of the present concealed state.

$$h_t=o_t\otimes \phi \left(c_t\right)$$

(32)

$$c_t=i_t\otimes c_t+f_t\otimes c_{t-1}$$

(33)

$$o_t=\sigma \left(W_o .\left[x_t.h_{t-1}\right]+b_o\right)$$

(34)

In the equation, $x_t$ represents the input, and $h_t$ signifies the network's output at time $t$. Meanwhile, $c_t$ represents the newly formed candidate value vectors for the tanh component at time $t$. The elements $f_t$, $i_t$, and $o_t$ pertain to the forget gate, input gates, and output gates, respectively. The matrices $W_f$, $W_i$, $W_0$, and $W_c$ correspond to the weights associated with the forget, input, output, and memory cells. Additionally, $b_f$, $b_i$, $b_o$, and $b_c$ denote the bias vectors linked to the forget, input, output, and memory cells, respectively. The function $\sigma$ symbolizes the sigmoid function, while $\phi$ is indicative of the Tanh function. To acquire the optimal parameters for augmenting the performance of the LSTM model, errors and gradients are calculated using the time backpropagation technique. The backpropagation in LSTM networks entails the computation of gradients over time, the modification of network parameters, and the control of information propagation through the cell state and gates. The LSTM's ability to capture long-term dependencies in sequential data makes it highly suitable for time-series prediction.

2.6. The Development of the Hybrid CEEMDAN-ARIMA-LSTM Model

In the development of the hybrid CEEMDAN-ARIMA-LSTM model, this work utilizes the complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN), autoregressive integrated moving average model (ARIMA), and long-term and short-term memory networks (LSTM). The CEEMDAN modal decomposition method is used to decompose the data at different frequencies to improve the timing characteristics reflected in the data and achieve the noise reduction effect, and the advantage of LSTM and ARIMA in time series data processing is exploited to enhance prediction accuracy. The specific procedures are listed below:

• Decompose the original data into Intrinsic Mode Functions (IMFs) and a residual component using the CEEMDAN technique. CEEMDAN decomposes the time series into numerous oscillatory mode components with varying frequencies, capturing both high-frequency and low-frequency components, thereby making it easier for subsequent models to capture distinct patterns.
• The ARIMA models are used to capture temporal dependencies and trends, as well as to analyse individually each of the retrieved IMFs and the residual component derived from CEEMDAN. Develop separate ARIMA models for the IMF and residual. The residual is then combined with the forecasts from the ARIMA models of each IMF to reconstruct the predicted series.
• From the ARIMA fitted results, the calculated residuals serve as inputs of the LSTM model. Train the LSTM model using the training set to generate 1-step ahead forecasts.
• The prediction result is derived by adding the predicted values of the high-frequency components using LSTM and the predicted value of the low-frequency components using ARIMA.
• Repeat steps (1) to (4) for $k=2, 3, \dots , K$ yields the final prediction.

2.7. Model Performance Criteria

The selection of the metric to assess the performance of the suggested model is of utmost importance. The selection of metrics has a significant impact on how the performance of model algorithms is evaluated and compared. The two most frequently used performance metrics are root mean square error (RMSE) and coefficient of determination $\left(R^2\right)$. RMSE measures the average difference between the actual and projected values, while $R^2$ quantifies the extent of correlation between them [42]. In addition, we utilize the concept of directional symmetry (DS) to evaluate a model's capacity to accurately predict the direction of changes in value. The metrics are mathematically defined as follows:

$$RMSE=\sqrt{{\displaystyle \frac{1}{n}}\sum _{i=1}^n{\left(y_i-{\widehat{y}}_i\right)}^2}$$

(35)

$$R^2=\sum _{i=1}^n{\left({\widehat{y}}_i-{\overline{y}}_i\right)}^2/\sum _{i=1}^n{\left(y_i-{\overline{y}}_i\right)}^2$$

(36)

$$DS\left(y,\widehat{y}\right)={\displaystyle \frac{100}{n-1}}\sum _{i=2}^n{d}_i\hspace{1em}\hspace{1em}\hspace{1em}{d}_i=\left\{\begin{array}{c}1, if \left(y_i-y_{i-1}\right)\left({\widehat{y}}_i-{\widehat{y}}_{i-1}\right)>0\\ 0, otherwise\hfill \end{array}\right.$$

(37)

where $n$ is the number of data points, $y_i$ is the observed value, ${\widehat{y}}_i$ is the predicted value and ${\overline{y}}_i$ is the mean of the observation data.

3. Results and Discussion

3.1. Rainfall Data Series and Trend Analysis

Figure 6 shows the monthly cumulative rainfall data measured at the Cape Town International Airport meteorological station in the Western Cape province of South Africa from 1995 to 2020. Typically, the monthly rainfall ranges from 0 to approximately 100 mm, with occasional instances where it exceeds 150 mm. The rainfall time series was utilised to calculate the Standardised Precipitation Index (SPI) in this research. Prior to calculating the SPI, it is crucial to evaluate the rainfall pattern across the research region. Thus, the section below will utilise innovative trend approaches to investigate rainfall trends.

The results of the ITA analysis of monthly total rainfall for the Cape Town International Airport are presented in Figure 7 and Table 2, estimated at a 95% confidence level. The 1:1 line of the ITA, with the 95% confidence limit area shaded light blue between the 95% confidence upper limit (blue line) and the 95% confidence lower limit (purple line). According to the ITA methodology, if the slope value falls within the lower and upper limits, there is no discernible trend. In contrast, if the slope exceeds the upper limit, an increasing trend is observed, whereas a decreasing trend is indicated if the slope falls below the lower limit. The rainfall variable results in Figure 7 revealed that most of the data points were situated below the 1:1 line, indicating a significant overall decreasing trend. This finding strongly aligned with the ITA trend indicator value (-3.215) (see Table 2). Additionally, Table 2 displays a slope of -0.083, which is below the ITA 95% confidence interval. This further confirms the presence of a decreasing trend, and the significance was established at the 5% level, as the calculated p-value was lower than the alpha level (0.05). However, it should be noted that the trend seems to recover for higher rainfall values.

The results of the sequential Mann-Kendall (SQ-MK) test statistics for the monthly rainfall dataset of Cape Town International. In this Figure monthly values of $u\left(t_i\right)$ and $u^{\prime }\left(t_i\right)$ for the period from 1995 to 2020. The major trend change detection point is reported in November, where the curves of $u\left(t_i\right)$ and $u^{\prime }\left(t_i\right) i$ntersect each other. This leads to a significant downwards trend that persists from 2010 onwards. The downward trend presented by SQ-MK seems to be consistent with the ITA results presented in Figure 7 and Table 2.

The MK and MMK techniques were also applied on the total rainfall time series (see Table 3). In Table 3, the results indicate that both the MK and the MMK methods reported a significant decreasing trend with negetive Z-score values of -3.374 and -4.077, respectively, and p-values of $1.74002\times {10}^{-04}$ and $4.5564\times {10}^{-05}$, respectively, which are lower than the alpha level (0.05). and in general, all these results indicate are consistent to those reported by the SQ-MK and ITA.

3.2. SPI Time Series and Forecusting Results

In this study, the SPI values for the 6-, 9-, and 12-month timescale were calculated using mothly total rainfall data measured at the Cape Town Internatinal Airport meteorological station in the Western Cape Region of South Africa (See Figure 6). Figure 9 illustrates the time series of SPI computed for timescales 6-months (SPI6), 9-months (SPI9), and 12-months (SPI12), respectively. In general, all SPIs (SPI6, SPI9, and SPI12) reported several episodes of moderate to high to excessive droughts in the study area. There is a broad peak of drought period that begins late 2009 and persist to 2012. SPI12 reported a persistant drought period that started in 2014, that led to day zero conditions in terms of water supply over Cape Town [43].

The complexity and fluctuation of the original SPI series in the model fitting and subsequent convergence of the model can be affected and thus limit the prediction accuracy. In response to the challenges, the nonlinear and nonstationary SPI series is preprocessed with the powerful method for conducting the decomposition of the series, namely CEEMDAN algorithm [44]. The decomposition of the time series is a fundamental component of the hybrid model proposed in this research for time series forecasting. As a result, CEEMDAN deconstructed the training dataset, yielding five IMF components and one trend time series. For an example of the decomposition of the time series, CEEMDAN was applied to the SPI6 time series, and the results are shown in Figure 10. While the CEEMDAN is able to present all IMFs, it is important to mention that IMF6 (which is assoiated with the Trend) seem to indicate that a downwards trend of the SPI6 time series, which could indicate that the domination of dry spells in recent over the study area.

This study compared the prediction performance of models listed in Table 4 before and after time series decomposition to determine whether the effort of decomposition improves the practicality of the model's prediction performance. Figure 11 presents a comparison of the prediction results of different models, along with the original time series for SPI6, SPI9, and SPI12, as well as a Taylor diagram. Overall, all models appear to closely mimic the original SPI time series across all time scales (see Figure 11). The Taylor diagram seem to indicate an improvement in prediction accuracy after applying CEEMDAN signal decomposition method, with the CEEMDAN-ARIMA-LSTM model outperforming other models in terms of prediction accuracy across all time scales. Table 4 evaluates the comparison of prediction performance values of different models using RMSE, $R^2$ and DS. As the time scale increases, the RMSE values of the models decrease, while the DS values generally increase (see Table 4). This indicates that the prediction accuracy of the models gradually improves with increasing time scale, peaking at the 12-month time scale. For instance, the LSTM model implemented at SPI6 had an RMSE of 0.234 and a $R^2$ of 0.897, while at SPI12, it had an RMSE of 0.058 and a $R^2$ of 0.984. In Table 4, it can be observed that at all-time scales, the RMSE values of the CEEMDAN-ARIMA and CEEMDAN-LSTM models were lower than those of the ARIMA and LSTM models, respectively, while the $R^2$ values were higher than those of the single models. This indicates higher prediction accuracy of the combined model, making it more suitable for predicting multiscale SPI. At each monthly time scale, the prediction accuracy of the ARIMA-LSTM combined model was significantly higher than that of the single model, with slightly higher accuracy for SPI6, SPI9, and SPI12. For example, in SPI6, the model had an RMSE of 0.186 and a $R^2$ of 0.931, while in SPI9, it had an RMSE of 0.077 and a $R^2$ of 0.983. In SPI12, the RMSE was 0.057 and the $R^2$was 0.985. It is evident that the prediction performance after the CEEMDAN decomposition is superior to that of the undecomposed models, suggesting that the SPI time series is better predicted after decomposition. Among these models, the CEEMDAN-ARIMA-LSTM model achieves the highest prediction accuracy with RMSE values ranging from 0.120 to 0.042, DS values ranging from 0.915 to 0.950, and values ranging from 0.970 to 0.995, significantly outperforming other models.

5. Conclusions

In the present study, the Mann-Kendall (MK) test and innovative trend analysis (ITA) method were utilized to analyze the trends in monthly rainfall total data at Cape Town International Airport, South Africa, for the period from 1995 to 2020. The results of the MK and MMK test indicated a significant decreasing trend in rainfall at a 95% significance level. These findings were further supported by the ITA results, with most data points falling below the 1:1 line, indicating a downward trend. Both the SQ-MK and ITA methods indicates a significant downwards trend of rainfall in the study area for most of the years. The recovery of the downwards trend seems to be observed during the period after the 2014-2016 strong El Niño years. Moreover, this study used rainfall data to compute the SPI and used this data train and forecast the SPI data at multiscale using ARIMA and LSTM. To achieve accurate SPI forecasts, the CEEMDAN signal processing algorithm was integrated with the ARIMA, LSTM, and ARIMA-LSTM hybrid model. A comprehensive comparative analysis of the prediction outcomes was conducted (see Figure 11 and Table 4). The results revealed that all models produced good performance in predicting SPI in all timescales. Performance steadily increased as the timescale increased presumably because of lower noise levels for higher SPI time scales. Throughout multiple timescales, the CEEMDAN combine hybrid model consistently outperformed the individual models. Specifically, for SPI12, the hybrid model exhibited RMSE values exceeding 0.97, highlighting CEEMDAN's ability to handle non-stationary and non-linear data characteristics, thereby improving predictability. Additionally, as depicted in Figure 11, the drought index predictions derived from the hybrid CEEMDAN-ARIMA-LSTM model are reported to have a superior prediction accuracy compared to other models across all timescales. It is evident that the combination of CEEMDAN with ARIMA, LSTM, and ARIMA-LSTM has the potential to significantly enhance the accuracy of meteorological drought forecasting.

The findings presented in this study seem to indicate that potential of hybrid models in predicting climate change patterns, which could be key for early warning systems. This could be useful for planning and policy decisions related to agricultural production and tourism management. Furthermore, analyzing climate parameters can help predict climate changes, especially drought and abnormal wet spells. And future investigations should explore the predictability of the standardized precipitation evapotranspiration index (SPEI) to assess the applicability of the combined model.