# Application of Extreme Learning Machine Algorithm for Drought Forecasting

Authors: Muhammad Ahmad Raza; Mohammed M. A. Almazah; Zulfiqar Ali; Ijaz Hussain; Fuad S. Al-Duais
Journal: Complexity
 (2022)

---

## Abstract

Drought is a complex and frequently occurring natural hazard in many parts of the world. Therefore, accurate drought forecasting is essential to mitigate its adverse impacts. This research has inferred the implication and the appropriateness of the extreme learning machine (ELM) algorithm for drought forecasting. For numerical evaluation, time series data of the Standardized Precipitating Temperature Index (SPTI) are used for nine meteorological stations located in various climatological zones of Pakistan. To assess the performance of ELM, this research includes parallel inferences of multilayer perceptron (MLP) and autoregressive integrated moving average (ARIMA) models. The performance of each model is assessed using root mean square error (RMSE), mean absolute error (MAE), mean absolute percent error (MAPE), Kling-Gupta efficiency (KGE), Willmott index (WI), and Karl Pearson’s correlation coefficient. Generally, graphical results illustrated an excellent performance of the ELM algorithm over MLP and ARIMA models. For training data of SPTI-1, ELM’s best performance has observed at Chitral station (RMSE = 0.374, KGE = 0.838, WI = 0.960, MAE = 0.272, MAPE = 259.59,R = 0.93). For SPTI-1 at Astore station, the numerical results are (RMSE = 0.688, KGE = 0.988, WI = 0.997, MAE = 0.798, MAPE = 247.35). The overall results indicate that the ELM outperformed by producing the smallest RMSE, MAE, and MAPE values and maximum values for KGE, WI, and correlation coefficient values at almost all the selected meteorological stations for (1, 3, 6, 9, and 12) month time scales. In summary, this research endorses the use of ELM for accurate drought forecasting.


---

## Body

## 1. Introduction

Drought is a recurrent natural climatic phenomenon that occurs virtually in most parts of the world. A drought is a recurrent event due to a lack of precipitation for an extended period of time in a particular region [1, 2]. Like other natural hazards, drought is steady and sometimes considered a creeping phenomenon as it is a gradually evolving natural hazard due to climatic fluctuations [3, 4]. Generally, the impacts of drought have effect on agriculture, livestock, ecological system, socio-economic, and energy sectors [5, 6]. Moreover, drought can be categorized as meteorological, hydrological, agricultural, and socio-economic drought. Therefore, it requires the investigator to consider individual opinions to define a specific type of drought [1, 2, 7, 8].Due to the complex nature of drought, it is difficult to monitor and assess its impact [1, 2, 6]. An accurate prediction of drought is considered difficult, especially its onset or end [9, 10]. Prolonged droughts adversely impact the economic agriculture and social sectors. These massive drought impacts are due to sudden and widespread climate changes [11]. Drought can lead to devastating economic effects, with worldwide losses of around $9 billion per annum; the US livestock industry faced a $400 million loss during a severe drought in 2002 [12]. A comprehensive early warning system for drought is necessary to reduce its devastating impact. However, a few studies are conducted to mitigate this stochastic natural hazard [13–15]. Recently, numerous drought indices have been developed to identify and monitor droughts and introduce mitigation policies [16–19]. Reference [20] proposed a new drought indicator, i.e., Normalized Ecosystem Drought Index (NEDI), to observe dryness conditions in the pattern of a transitional ecosystem. It is expected that dryness conditions can be quantified better by using NEDI.Numerically expressed drought indices are more understandable than natural rainfall data [1, 2, 21]. Drought indices can be a valuable tool to detect the initiation and termination of drought levels necessary for recovery planning, mitigation, and decision-making [22, 23]. Drought indices aim to quantify how drought conditions evolve and classify the severity of drought events. These indices made easy droughts modeling using stochastic time series, neural network algorithms, and water balance models. The most commonly used drought indices are Palmer Drought Severity Index (PDSI) [24], Surface Water Supply Index (SWSI) [25], Standardized Precipitation Index (SPI) [26], Effective Drought Index (EDI), and Standardized Precipitation and Evapotranspiration Index (SPEI) [27]. Ali et al. [28] proposed a multiscalar drought index named as Standardized Precipitation Temperature Index (SPTI). These drought indices were calculated using different meteorological variables [29]. Different drought indices were used to characterize, estimate, and forecast drought conditions.The current long-range drought forecasts have minimal reliability [30]. Existing conventional stochastic models are inadequate for accurate drought predictions [31]. The recently developed machine learning (ML) models have extensive application in climatology including Naïve Bayes classifier, Bayesian networks [32], support vector machine (SVM), wavelet gene expression programming [33], maximum entropy, and artificial neural networks (ANNs). Results of several studies affirmed that the ML models perform comparatively better than conventional stochastic and dynamic models for drought estimation [34, 35]. ANN models act like a human brain and can be classified according to their neuron structure, number of hidden layers, and activation functions. Many researchers have successfully applied MLP neural networks for drought estimation and forecasting [36, 37]. The MLP is capable of accurately forecasting soil temperature in semi-humid and arid regions [38]. Aghelpour et al. [39] improved agriculture drought modeling by coupling the dragonfly optimization algorithm with SVM. Furthermore, [40] efficiently modeled RDI using hybrid support vector regression (SVR) coupled with firefly algorithm (FA), whale optimization algorithm (WOA), and wavelet analysis (WA). The results proved that hybrid and coupled SVR techniques improved drought forecasting. Although ML models have an outstanding reputation in estimation, prediction, and forecasting, many have slow computing times [41]. Among the class of ANN algorithms, ELM is being widely used in various fields and has gained fame in climatology and engineering [42–47]. Mouatadid and Adamowski [48] efficiently forecasted urban water demand for Montreal (Canada) using ELM.This research aims to infer the implication and the appropriateness of the extreme learning machine (ELM) algorithm for drought forecasting. In previous research, ELM has been implemented in different disciplines, including classification [49], regression [50, 51], clustering [52, 53], feature selection [54], pattern recognition [55], image processing [56] estimating sediment transport [57], and drought forecasting [30, 58, 59].The ELM has significantly faster learning, improved generalization performance, minimum human intervention, and accurate forecasting performance [60].


## 2. Materials and Methods

### 2.1. Data and Study Area

The application of this research is based on nine meteorological stations scattered around Pakistan. The topographic map of the study region and distribution of selected meteorological stations is shown in Figure1. The study area is situated in the southeastern part of Asia and lies between 23.8° to 37°N latitude and 60.9° to 75.37°E longitude. The region is classified into clusters comprising different meteorological stations with diverse spatial characteristics [61]. Hence, selecting these meteorological stations aims to cover the maximum climatic variability. In addition, the study area encompasses five major river basins, Ravi, Chenab, Sutlej, Jhelum, and the Indus River. These rivers are the backbone of the country’s agriculture industry and hydropower projects.Figure 1
Topographic map of the study area and geographical distribution of meteorological stations.For this research, time series data of the monthly precipitation and minimum and maximum air temperatures were collected from the Karachi Data Processing Center (KDPC) through the Pakistan Meteorological Department (PMD). The length of the data ranges from January 1951 to December 2016. The full-length data were split into two parts. January 1951 to December 2013 is considered the training data set, and the remaining three years, January 2014 to December 2016, is considered the test data set. The climatological forecast needs more accuracy for future hazard mitigation because long-range climatological forecasts compromise accuracy. Here, the errors and irregularities were detected and removed by the KDPC itself. Additionally, missing data were adjusted by generating values using cumulative distributions over lead periods.


### 2.2. Standardized Precipitation Temperature Index (SPTI)

The Standardized Drought Indices (SDIs) have extensive applications for drought monitoring. SDIs are standardized and spatially invariant tools for monitoring and assessing drought characteristics. In the literature, various authors have offered numerous methods for SDIs. Example includes the Standardized Precipitation Index (SPI) [26], Standardized Precipitation Evapotranspiration Index (SPEI) [27], and Standardized Precipitation Temperature Index (SPTI) [28]. Precipitation and temperature are two essential climatology indicators, revealing the vital dynamics of climate and hydrology. Therefore, a standardized drought index based on these two meteorological variables is more beneficial for drought monitoring and forecasting. Therefore, the SPTI has been chosen as SDI for monitoring and forecasting drought. The mathematical calculation of SPTI is quite similar to SPI; more detailed discussion can be accessed in [41]. SPTI is a multiscalar drought index and can be calculated for different time scales (1–48). Positive and negative values of the index indicate drought and wet conditions. These drought conditions are classified in Table 1 [62, 63].Table 1
Drought classification of the SPTI index.
SPTI valuesClassesCumulative probabilities≤−2Extreme drought (ED)0.000–0.023−1.9 to −1.5Severe drought (SD)0.023–0.067−1.49 to −1Moderate drought (MD)0.067–0.159−0.99 to 0.99Near normal (NN)0.159–0.8411 to 1.49Moderate wet (MW)0.841–0.9331.5 to 1.9Severe wet (SW)0.933–0.977≥2Extreme wet (EW)0.977–1SPTI is a modified form of the De-Martonne Aridity Index (DAI) (de Martonne, 1926). The mathematical properties of SPTI are utterly similar to SPI, an extensively used index for drought prediction in many parts of the world. For SPTI, we need to calculate DAI based on the monthly total precipitation and average monthly temperature. The next step is to fit an appropriate distribution to calculate a cumulative probability for standardization. However, many researchers used Gamma distribution for standardization. The index values are subjected to fitted distribution, and none of the single distribution can be appropriate for all the stations and for various time scales. Therefore, the 32 candidate distributions have been fitted on DAI at different lead time scales. The Bayesian Information Criterion (BIC) has been used as a threshold to assess the appropriateness of a distribution.


### 2.3. Candidate Algorithms

An artificial neural network (ANN) is a computational paradigm. It is a data-driven technique in which information goes through a biological structure of neurons with multiple layers introduced in the 1950s. It did not impose any constraints on input variables to train the model like other stochastic models. These algorithms are brilliant and learn from existing relationships among the observations of input and auxiliary variables. ANN can manage high-dimensional and high-frequency complex datasets [58]. ANN algorithms have broad applications in mathematics, engineering, medicine, economics, neurology, and hydrology [64–67]. Kuligowski and Barros [68] claimed that weather prediction could be improved using ANN algorithms. This class of algorithms can be helpful in the field of climatology to forecast natural hazards like drought. Multilayer perceptron (MLP) is considered one of the useful and fully connected feedforward artificial neural networks. It usually consists of three layers of multiple nodes, including an input layer, multiple hidden layers (usually two hidden layers with multiple hidden nodes) with a nonlinear activation function, and an output layer. The neuron structure and estimation accuracy of MLP make it prominent among the other ANN algorithms. Error backpropagation is one of the supervised learning techniques used to train MLP. Another stochastic algorithm used for drought prediction is the ARIMA process [69].


#### 2.3.1. Seasonal Autoregressive and Integrated Moving Average Model (SARIMA)

Yule [70] pioneered to introduce autoregressive (AR) models that the time series being analyzed is a linear function of its previous lag values. Slutzky [71] modeled time series as a function of past residual terms named as moving average model (MA). Wold [72] merged both AR and MA specifications and introduced a new generalized form of ARMA specifications used to model all stationary time series by choosing the appropriate order of “AR” and “MA” terms into the model. Time series data generally have trends (non-stationary). Non-stationary time series can be modeled by appropriate differencing the series into stationary. The series that is transformed from non-stationary to stationary by differencing is known as integrated series. The ARIMA has a systematic way of identification, estimation, and diagnostic checking approach to reach an appropriate model. Many hydrologic and meteorological time series data have inherited seasonal components [73]. These kinds of data can be efficiently modeled with the seasonal ARIMA model, which requires only a few parameters to be estimated [74]. The seasonal ARIMA model is described as ARIMA(p, d, q) (P, D, Q)s, where (p, d, q) is the non-seasonal component of the ARIMA specifications, while (P,D,Q)s is the seasonal component of the ARIMA model. The general seasonal ARIMA specifications are as follows:(1)ϕpBΦPBs∇d∇sD′Dt=θqBΘQBsαt.Here, “p” is the order of non-seasonal autoregressive terms, “q” is the no of non-seasonal MA terms to be included in the model. Similarly, “P” is the seasonal autoregressive terms, and “Q” is the number of seasonal MA terms, ∇d is the difference operator of non-seasonal series with “d” levels to make the series stationary, ∇D′ is the difference operator of seasonal series with D′ no of differencing to get integrated stationary series, where “s” is the length of the season. The mathematical details of ARIMA specifications can be observed in [75]. The development of ARIMA specifications included identification, estimation, and diagnostic checking. By following these steps, a parsimonious ARIMA specification can be selected for estimation and forecasting a time series.


#### 2.3.2. Extreme Learning Machine (ELM)

The extreme learning machine (ELM) is a modern single hidden layer feedforward neural network (SLFN) algorithm proposed by [76]. The proposed novel machine learning algorithm (ELM) operates identically to feedforward back-propagation ANN (FFBP-ANN) and least-squares support vector regression LSSVR models. It has shown its candidacy among the ANN algorithms to solve complex linear and nonlinear regression problems. It contains a single hidden layer of multiple hidden nods. However, most of the ANN-based methods have specific limitations such as slow computation, learning epochs, larger biases, and tuning parameters (weights). To overcome such weaknesses and the frailty of ANN methods, a state-of-the-art algorithm known as extreme learning machine (ELM) gained fame in the class of ANN algorithms [77].Studies have revealed that even with randomly generated weights of hidden nodes, ELM can attain the universal approximation feature of SLFNs [46, 47, 78]. In the proposed method, the input weights are assigned randomly, and the output weights can be solved uniquely by the least-squares method of generalized inverse function [76]. If hidden node input weights and biases are chosen randomly, SLNFs can be considered a linear system. The output weights are determined analytically through the generalized inverse operation of the hidden layer output matrices because these weights connect the hidden layer to the output layer of the linear system. The ELM can solve regression problems with shorter simulation times than FFBP-ANN and LSSVR algorithms and makes ELM a thousand times faster [59, 79–81]. It contains common properties of high generalization performance. The topological structure of the ELM algorithm is given in Figure 2, where three layers of neurons are used to develop the architecture.Figure 2
The topological structure of the extreme learning machine network. The input layer is the lag values of SPTI, the hidden layer contains randomly assigned hidden nodes, and the output layer generates predicted values of SPTI.Input layer where the input variables are introduced, the single hidden layer contain variable number of neurons where data are processed and analyzed and the output layer produce desired results through their activation function. The activation function used in the ELM algorithm is a sigmoid functiongx=1/1+e−x, and while training ELM, most of the time is utilized while calculating the Moore-Penrose generalized inverse H† of the hidden layer.Huang et al. [76] claimed that the maximal margin property of SVM and the minimal norm of weight theory of ELM is consistent. ELM and SVM perform equally well for standard optimization. For M random distinct samples yi,ti, where yi=yi1,yi2,…,yimT∈Rm and ti=ti1,ti2,…tilT∈Rl, standard ELM with N hidden nodes and with activation function gy are mathematically modeled as(2)∑i=1Ν⌢γigiyj=∑i=1Ν⌢γigωi.yj+bi=ojj=1,2,3,…,M,where ωi=ωi1,ωi2,…,ωimT is the weight vector connecting the ith hidden node and input nodes. γi=γi1,γi2,…,γilT is the weight vector connecting the ith hidden node and output nodes, and bi is the threshold of ith hidden node. ωi.yj denote the inner product of ωi and yj.(3)Hγ=T,(4)Hω1,ω2,…ωM⌢,b1,b2,…,bM⌢,y1,y2,…,yM=gω1.y1+b1⋯gωM⌢.y1+bM⌢⋮⋯⋮gω1.yM+b1⋯gωM⌢.yM+bM⌢M⌢×M,(5)γ=γ1Tγ2T⋮γM⌢TM⌢×l,andT=t1Tt2T⋮tMTM×l.ELM attains optimal generalization performance as long as the chosen number of hidden nodes is sufficiently high. In our simulation through ELM with sigmoid activation function, the number of hidden nodes is selected automatically to attain optimal prediction and forecast performance.


### 2.4. Model Evaluation Metrics

A model performance assessment needs calibration of an existed link between observed and predicted hydrological patterns. The fundamental performance assessment method is through a visual inspection of empirical and predicted or forecasted time series. For the quantitative evaluation of algorithms, around 20 performance metrics select a hydrological model [82]. It has been observed that the choice of an appropriate model significantly changed if precision-based metrics were used instead of error-based metrics [83]. Numerous accuracy measurement criteria were developed, but each tool has inherited pros and cons, and none of the metrics is universally accepted and can be used as a threshold [84]. In this study, some error-based performance metrics have been used for cross-validation, including root mean square error (RMSE), mean absolute error (MAE), and mean absolute percentage error (MAPE). The Kling-Gupta efficiency (KGE) and Willmott Index of agreement (WI) are also better ways to assess the performance of stochastic, machine learning, and hydrologic models [85, 86]. Another way to assess an algorithm’s prediction performance is to calculate the simple correlation coefficient between observed and predicted values of the input variable (SPTI) as a closeness measure. Similar performance metrics are used to assess the forecast ability of candidate algorithms.The RMSE is the deviation of estimated or predicted values “D^” from actual or observed values “D” of drought indices, computed for “T” different predictions given in(6)RMSE=1h∑t=TT+hDt+1−D^t+1|T2.Since RMSE is positively affected by outliers, therefore we need some robust measures toward extreme values. Mean absolute error (MAE) is less influenced by extreme values than RMSE [87]. Equation (7) describes the mathematical structure of the MAE.(7)MAE=1h∑t=TT+hDt+1−D^t+1|T.Another accuracy measure is the mean absolute percentage error (MAPE), a unit-free tool to assess an algorithm’s prediction and forecast ability. Unlike other performance metrics, MAPE is a scaled independent metric. These performance metrics or accuracy measures are extensively being used in the field of climatology. The mathematical form of MAPE is given in(8)MAPE=100h∑t=TT+hDt+1−D^t+1|TDt+1.Kling-Gupta efficiency index was developed to assess the performance model by comparing estimated and observed time series data [88].(9)KGE=1−r−12+α−12+β−12.Here, “r,” α, and β in the KGE index illustrate the correlation coefficient, standard deviation ratio, and average ratio of observed and predicted values of SPTI, respectively.Willmott [89] proposes an index named Willmott Index of agreement (WI) as a standardized measure of the degree of model prediction error.(10)WI=1−∑t=TT+hDt+1−D^t+1|T2∑t=TT+hDt+1−D¯t+1+D^t+1−D¯t+12,0≤WI≤1.A model with minimum values of RMSE, MAE, MAPE, the maximum value of KGE index, and the value of WI close to “1” will be selected and proposed as an adequate algorithm for the estimation of existing drought phenomena and forecasting future drought episodes.


## 2.1. Data and Study Area

The application of this research is based on nine meteorological stations scattered around Pakistan. The topographic map of the study region and distribution of selected meteorological stations is shown in Figure1. The study area is situated in the southeastern part of Asia and lies between 23.8° to 37°N latitude and 60.9° to 75.37°E longitude. The region is classified into clusters comprising different meteorological stations with diverse spatial characteristics [61]. Hence, selecting these meteorological stations aims to cover the maximum climatic variability. In addition, the study area encompasses five major river basins, Ravi, Chenab, Sutlej, Jhelum, and the Indus River. These rivers are the backbone of the country’s agriculture industry and hydropower projects.Figure 1
Topographic map of the study area and geographical distribution of meteorological stations.For this research, time series data of the monthly precipitation and minimum and maximum air temperatures were collected from the Karachi Data Processing Center (KDPC) through the Pakistan Meteorological Department (PMD). The length of the data ranges from January 1951 to December 2016. The full-length data were split into two parts. January 1951 to December 2013 is considered the training data set, and the remaining three years, January 2014 to December 2016, is considered the test data set. The climatological forecast needs more accuracy for future hazard mitigation because long-range climatological forecasts compromise accuracy. Here, the errors and irregularities were detected and removed by the KDPC itself. Additionally, missing data were adjusted by generating values using cumulative distributions over lead periods.


## 2.2. Standardized Precipitation Temperature Index (SPTI)

The Standardized Drought Indices (SDIs) have extensive applications for drought monitoring. SDIs are standardized and spatially invariant tools for monitoring and assessing drought characteristics. In the literature, various authors have offered numerous methods for SDIs. Example includes the Standardized Precipitation Index (SPI) [26], Standardized Precipitation Evapotranspiration Index (SPEI) [27], and Standardized Precipitation Temperature Index (SPTI) [28]. Precipitation and temperature are two essential climatology indicators, revealing the vital dynamics of climate and hydrology. Therefore, a standardized drought index based on these two meteorological variables is more beneficial for drought monitoring and forecasting. Therefore, the SPTI has been chosen as SDI for monitoring and forecasting drought. The mathematical calculation of SPTI is quite similar to SPI; more detailed discussion can be accessed in [41]. SPTI is a multiscalar drought index and can be calculated for different time scales (1–48). Positive and negative values of the index indicate drought and wet conditions. These drought conditions are classified in Table 1 [62, 63].Table 1
Drought classification of the SPTI index.
SPTI valuesClassesCumulative probabilities≤−2Extreme drought (ED)0.000–0.023−1.9 to −1.5Severe drought (SD)0.023–0.067−1.49 to −1Moderate drought (MD)0.067–0.159−0.99 to 0.99Near normal (NN)0.159–0.8411 to 1.49Moderate wet (MW)0.841–0.9331.5 to 1.9Severe wet (SW)0.933–0.977≥2Extreme wet (EW)0.977–1SPTI is a modified form of the De-Martonne Aridity Index (DAI) (de Martonne, 1926). The mathematical properties of SPTI are utterly similar to SPI, an extensively used index for drought prediction in many parts of the world. For SPTI, we need to calculate DAI based on the monthly total precipitation and average monthly temperature. The next step is to fit an appropriate distribution to calculate a cumulative probability for standardization. However, many researchers used Gamma distribution for standardization. The index values are subjected to fitted distribution, and none of the single distribution can be appropriate for all the stations and for various time scales. Therefore, the 32 candidate distributions have been fitted on DAI at different lead time scales. The Bayesian Information Criterion (BIC) has been used as a threshold to assess the appropriateness of a distribution.


## 2.3. Candidate Algorithms

An artificial neural network (ANN) is a computational paradigm. It is a data-driven technique in which information goes through a biological structure of neurons with multiple layers introduced in the 1950s. It did not impose any constraints on input variables to train the model like other stochastic models. These algorithms are brilliant and learn from existing relationships among the observations of input and auxiliary variables. ANN can manage high-dimensional and high-frequency complex datasets [58]. ANN algorithms have broad applications in mathematics, engineering, medicine, economics, neurology, and hydrology [64–67]. Kuligowski and Barros [68] claimed that weather prediction could be improved using ANN algorithms. This class of algorithms can be helpful in the field of climatology to forecast natural hazards like drought. Multilayer perceptron (MLP) is considered one of the useful and fully connected feedforward artificial neural networks. It usually consists of three layers of multiple nodes, including an input layer, multiple hidden layers (usually two hidden layers with multiple hidden nodes) with a nonlinear activation function, and an output layer. The neuron structure and estimation accuracy of MLP make it prominent among the other ANN algorithms. Error backpropagation is one of the supervised learning techniques used to train MLP. Another stochastic algorithm used for drought prediction is the ARIMA process [69].


### 2.3.1. Seasonal Autoregressive and Integrated Moving Average Model (SARIMA)

Yule [70] pioneered to introduce autoregressive (AR) models that the time series being analyzed is a linear function of its previous lag values. Slutzky [71] modeled time series as a function of past residual terms named as moving average model (MA). Wold [72] merged both AR and MA specifications and introduced a new generalized form of ARMA specifications used to model all stationary time series by choosing the appropriate order of “AR” and “MA” terms into the model. Time series data generally have trends (non-stationary). Non-stationary time series can be modeled by appropriate differencing the series into stationary. The series that is transformed from non-stationary to stationary by differencing is known as integrated series. The ARIMA has a systematic way of identification, estimation, and diagnostic checking approach to reach an appropriate model. Many hydrologic and meteorological time series data have inherited seasonal components [73]. These kinds of data can be efficiently modeled with the seasonal ARIMA model, which requires only a few parameters to be estimated [74]. The seasonal ARIMA model is described as ARIMA(p, d, q) (P, D, Q)s, where (p, d, q) is the non-seasonal component of the ARIMA specifications, while (P,D,Q)s is the seasonal component of the ARIMA model. The general seasonal ARIMA specifications are as follows:(1)ϕpBΦPBs∇d∇sD′Dt=θqBΘQBsαt.Here, “p” is the order of non-seasonal autoregressive terms, “q” is the no of non-seasonal MA terms to be included in the model. Similarly, “P” is the seasonal autoregressive terms, and “Q” is the number of seasonal MA terms, ∇d is the difference operator of non-seasonal series with “d” levels to make the series stationary, ∇D′ is the difference operator of seasonal series with D′ no of differencing to get integrated stationary series, where “s” is the length of the season. The mathematical details of ARIMA specifications can be observed in [75]. The development of ARIMA specifications included identification, estimation, and diagnostic checking. By following these steps, a parsimonious ARIMA specification can be selected for estimation and forecasting a time series.


### 2.3.2. Extreme Learning Machine (ELM)

The extreme learning machine (ELM) is a modern single hidden layer feedforward neural network (SLFN) algorithm proposed by [76]. The proposed novel machine learning algorithm (ELM) operates identically to feedforward back-propagation ANN (FFBP-ANN) and least-squares support vector regression LSSVR models. It has shown its candidacy among the ANN algorithms to solve complex linear and nonlinear regression problems. It contains a single hidden layer of multiple hidden nods. However, most of the ANN-based methods have specific limitations such as slow computation, learning epochs, larger biases, and tuning parameters (weights). To overcome such weaknesses and the frailty of ANN methods, a state-of-the-art algorithm known as extreme learning machine (ELM) gained fame in the class of ANN algorithms [77].Studies have revealed that even with randomly generated weights of hidden nodes, ELM can attain the universal approximation feature of SLFNs [46, 47, 78]. In the proposed method, the input weights are assigned randomly, and the output weights can be solved uniquely by the least-squares method of generalized inverse function [76]. If hidden node input weights and biases are chosen randomly, SLNFs can be considered a linear system. The output weights are determined analytically through the generalized inverse operation of the hidden layer output matrices because these weights connect the hidden layer to the output layer of the linear system. The ELM can solve regression problems with shorter simulation times than FFBP-ANN and LSSVR algorithms and makes ELM a thousand times faster [59, 79–81]. It contains common properties of high generalization performance. The topological structure of the ELM algorithm is given in Figure 2, where three layers of neurons are used to develop the architecture.Figure 2
The topological structure of the extreme learning machine network. The input layer is the lag values of SPTI, the hidden layer contains randomly assigned hidden nodes, and the output layer generates predicted values of SPTI.Input layer where the input variables are introduced, the single hidden layer contain variable number of neurons where data are processed and analyzed and the output layer produce desired results through their activation function. The activation function used in the ELM algorithm is a sigmoid functiongx=1/1+e−x, and while training ELM, most of the time is utilized while calculating the Moore-Penrose generalized inverse H† of the hidden layer.Huang et al. [76] claimed that the maximal margin property of SVM and the minimal norm of weight theory of ELM is consistent. ELM and SVM perform equally well for standard optimization. For M random distinct samples yi,ti, where yi=yi1,yi2,…,yimT∈Rm and ti=ti1,ti2,…tilT∈Rl, standard ELM with N hidden nodes and with activation function gy are mathematically modeled as(2)∑i=1Ν⌢γigiyj=∑i=1Ν⌢γigωi.yj+bi=ojj=1,2,3,…,M,where ωi=ωi1,ωi2,…,ωimT is the weight vector connecting the ith hidden node and input nodes. γi=γi1,γi2,…,γilT is the weight vector connecting the ith hidden node and output nodes, and bi is the threshold of ith hidden node. ωi.yj denote the inner product of ωi and yj.(3)Hγ=T,(4)Hω1,ω2,…ωM⌢,b1,b2,…,bM⌢,y1,y2,…,yM=gω1.y1+b1⋯gωM⌢.y1+bM⌢⋮⋯⋮gω1.yM+b1⋯gωM⌢.yM+bM⌢M⌢×M,(5)γ=γ1Tγ2T⋮γM⌢TM⌢×l,andT=t1Tt2T⋮tMTM×l.ELM attains optimal generalization performance as long as the chosen number of hidden nodes is sufficiently high. In our simulation through ELM with sigmoid activation function, the number of hidden nodes is selected automatically to attain optimal prediction and forecast performance.


## 2.3.1. Seasonal Autoregressive and Integrated Moving Average Model (SARIMA)

Yule [70] pioneered to introduce autoregressive (AR) models that the time series being analyzed is a linear function of its previous lag values. Slutzky [71] modeled time series as a function of past residual terms named as moving average model (MA). Wold [72] merged both AR and MA specifications and introduced a new generalized form of ARMA specifications used to model all stationary time series by choosing the appropriate order of “AR” and “MA” terms into the model. Time series data generally have trends (non-stationary). Non-stationary time series can be modeled by appropriate differencing the series into stationary. The series that is transformed from non-stationary to stationary by differencing is known as integrated series. The ARIMA has a systematic way of identification, estimation, and diagnostic checking approach to reach an appropriate model. Many hydrologic and meteorological time series data have inherited seasonal components [73]. These kinds of data can be efficiently modeled with the seasonal ARIMA model, which requires only a few parameters to be estimated [74]. The seasonal ARIMA model is described as ARIMA(p, d, q) (P, D, Q)s, where (p, d, q) is the non-seasonal component of the ARIMA specifications, while (P,D,Q)s is the seasonal component of the ARIMA model. The general seasonal ARIMA specifications are as follows:(1)ϕpBΦPBs∇d∇sD′Dt=θqBΘQBsαt.Here, “p” is the order of non-seasonal autoregressive terms, “q” is the no of non-seasonal MA terms to be included in the model. Similarly, “P” is the seasonal autoregressive terms, and “Q” is the number of seasonal MA terms, ∇d is the difference operator of non-seasonal series with “d” levels to make the series stationary, ∇D′ is the difference operator of seasonal series with D′ no of differencing to get integrated stationary series, where “s” is the length of the season. The mathematical details of ARIMA specifications can be observed in [75]. The development of ARIMA specifications included identification, estimation, and diagnostic checking. By following these steps, a parsimonious ARIMA specification can be selected for estimation and forecasting a time series.


## 2.3.2. Extreme Learning Machine (ELM)

The extreme learning machine (ELM) is a modern single hidden layer feedforward neural network (SLFN) algorithm proposed by [76]. The proposed novel machine learning algorithm (ELM) operates identically to feedforward back-propagation ANN (FFBP-ANN) and least-squares support vector regression LSSVR models. It has shown its candidacy among the ANN algorithms to solve complex linear and nonlinear regression problems. It contains a single hidden layer of multiple hidden nods. However, most of the ANN-based methods have specific limitations such as slow computation, learning epochs, larger biases, and tuning parameters (weights). To overcome such weaknesses and the frailty of ANN methods, a state-of-the-art algorithm known as extreme learning machine (ELM) gained fame in the class of ANN algorithms [77].Studies have revealed that even with randomly generated weights of hidden nodes, ELM can attain the universal approximation feature of SLFNs [46, 47, 78]. In the proposed method, the input weights are assigned randomly, and the output weights can be solved uniquely by the least-squares method of generalized inverse function [76]. If hidden node input weights and biases are chosen randomly, SLNFs can be considered a linear system. The output weights are determined analytically through the generalized inverse operation of the hidden layer output matrices because these weights connect the hidden layer to the output layer of the linear system. The ELM can solve regression problems with shorter simulation times than FFBP-ANN and LSSVR algorithms and makes ELM a thousand times faster [59, 79–81]. It contains common properties of high generalization performance. The topological structure of the ELM algorithm is given in Figure 2, where three layers of neurons are used to develop the architecture.Figure 2
The topological structure of the extreme learning machine network. The input layer is the lag values of SPTI, the hidden layer contains randomly assigned hidden nodes, and the output layer generates predicted values of SPTI.Input layer where the input variables are introduced, the single hidden layer contain variable number of neurons where data are processed and analyzed and the output layer produce desired results through their activation function. The activation function used in the ELM algorithm is a sigmoid functiongx=1/1+e−x, and while training ELM, most of the time is utilized while calculating the Moore-Penrose generalized inverse H† of the hidden layer.Huang et al. [76] claimed that the maximal margin property of SVM and the minimal norm of weight theory of ELM is consistent. ELM and SVM perform equally well for standard optimization. For M random distinct samples yi,ti, where yi=yi1,yi2,…,yimT∈Rm and ti=ti1,ti2,…tilT∈Rl, standard ELM with N hidden nodes and with activation function gy are mathematically modeled as(2)∑i=1Ν⌢γigiyj=∑i=1Ν⌢γigωi.yj+bi=ojj=1,2,3,…,M,where ωi=ωi1,ωi2,…,ωimT is the weight vector connecting the ith hidden node and input nodes. γi=γi1,γi2,…,γilT is the weight vector connecting the ith hidden node and output nodes, and bi is the threshold of ith hidden node. ωi.yj denote the inner product of ωi and yj.(3)Hγ=T,(4)Hω1,ω2,…ωM⌢,b1,b2,…,bM⌢,y1,y2,…,yM=gω1.y1+b1⋯gωM⌢.y1+bM⌢⋮⋯⋮gω1.yM+b1⋯gωM⌢.yM+bM⌢M⌢×M,(5)γ=γ1Tγ2T⋮γM⌢TM⌢×l,andT=t1Tt2T⋮tMTM×l.ELM attains optimal generalization performance as long as the chosen number of hidden nodes is sufficiently high. In our simulation through ELM with sigmoid activation function, the number of hidden nodes is selected automatically to attain optimal prediction and forecast performance.


## 2.4. Model Evaluation Metrics

A model performance assessment needs calibration of an existed link between observed and predicted hydrological patterns. The fundamental performance assessment method is through a visual inspection of empirical and predicted or forecasted time series. For the quantitative evaluation of algorithms, around 20 performance metrics select a hydrological model [82]. It has been observed that the choice of an appropriate model significantly changed if precision-based metrics were used instead of error-based metrics [83]. Numerous accuracy measurement criteria were developed, but each tool has inherited pros and cons, and none of the metrics is universally accepted and can be used as a threshold [84]. In this study, some error-based performance metrics have been used for cross-validation, including root mean square error (RMSE), mean absolute error (MAE), and mean absolute percentage error (MAPE). The Kling-Gupta efficiency (KGE) and Willmott Index of agreement (WI) are also better ways to assess the performance of stochastic, machine learning, and hydrologic models [85, 86]. Another way to assess an algorithm’s prediction performance is to calculate the simple correlation coefficient between observed and predicted values of the input variable (SPTI) as a closeness measure. Similar performance metrics are used to assess the forecast ability of candidate algorithms.The RMSE is the deviation of estimated or predicted values “D^” from actual or observed values “D” of drought indices, computed for “T” different predictions given in(6)RMSE=1h∑t=TT+hDt+1−D^t+1|T2.Since RMSE is positively affected by outliers, therefore we need some robust measures toward extreme values. Mean absolute error (MAE) is less influenced by extreme values than RMSE [87]. Equation (7) describes the mathematical structure of the MAE.(7)MAE=1h∑t=TT+hDt+1−D^t+1|T.Another accuracy measure is the mean absolute percentage error (MAPE), a unit-free tool to assess an algorithm’s prediction and forecast ability. Unlike other performance metrics, MAPE is a scaled independent metric. These performance metrics or accuracy measures are extensively being used in the field of climatology. The mathematical form of MAPE is given in(8)MAPE=100h∑t=TT+hDt+1−D^t+1|TDt+1.Kling-Gupta efficiency index was developed to assess the performance model by comparing estimated and observed time series data [88].(9)KGE=1−r−12+α−12+β−12.Here, “r,” α, and β in the KGE index illustrate the correlation coefficient, standard deviation ratio, and average ratio of observed and predicted values of SPTI, respectively.Willmott [89] proposes an index named Willmott Index of agreement (WI) as a standardized measure of the degree of model prediction error.(10)WI=1−∑t=TT+hDt+1−D^t+1|T2∑t=TT+hDt+1−D¯t+1+D^t+1−D¯t+12,0≤WI≤1.A model with minimum values of RMSE, MAE, MAPE, the maximum value of KGE index, and the value of WI close to “1” will be selected and proposed as an adequate algorithm for the estimation of existing drought phenomena and forecasting future drought episodes.


## 3. Results

The descriptive statistics of meteorological and climatic variables are briefly detailed by using five-number summary statistics. The numerical results related to Minimum (Min.), first quartile (Q1), Median, Mean, and third quartile (Q3) are expressed in Table 2. These results indicated that the annual and seasonal meteorological characteristics of the selected meteorological stations are quite diverse.Table 2
Summary statistics of rainfall (in mm) and spatial characteristics of selected stations.
StationMinQ1MedianMeanQ3MaxLatitudeLongitudeTime rangeAstore010.9726.3540.6454.52248.735.3570°E74.8624°N1954–2016Chhor00018.8911.4738125.5114°E69.7823°N1951–2016Chitral03.7220.238.5158.48281.435.7699°E71.7741°N1965–2016Kalat001.3515.3516.23546.029.0523°E66.5879°N1951–2016Kohat012.513347.9865.7419.633.5889°E71.4429°N1954–2016Mianwali03.9522.143.660.1253032.5839°E71.5370°N1960–2016Multan00.014.116.8519.2231.230.1575°E71.5247°N1951–2016Muzaffarabad046.3397.85125.65178.0772134.3551°E73.4769°N1955–2016Sialkot08.2833.0580.69296.3917.632.4945°E74.5229°N1951–2016Muzaffarabad has the highest mean monthly precipitation (125.65 mm), and the lowest mean monthly precipitation recorded was (15.35 mm) at Kalat. Sialkot has the highest maximum rainfall in a month (917.6 mm), and Chhor has observed the lowest maximum rainfall in a month (11.47 mm), while minimum rainfall at all selected locations was zero (0 mm). The precipitation source at these stations varies, such as heavy rainfall occurring at some stations in the monsoon season (June-Sep). However, precipitation exponentially declined after September and lasted till December until the western depression started in winter. At the same time, western depression causes rainfall in the winter season (Dec-Mar). The above statistics exhibit the dry and wet season cycles at a few stations. Temperature is another climatic variable used to calculate SPTI, so similar descriptive statistics for minimum and maximum temperature are expressed in Table3. Results show that the total monthly minimum and maximum temperatures are highly apparent and distinguishable for all the stations.Table 3
Summary statistics of minimum and maximum temperatures (Celsius) at selected meteorological stations.
StationMinimum temperatureMaximum temperatureElevation (m)MinQ1MedianMeanQ3MaxMinQ1MedianMeanQ3MaxAstore−12.1−2.54.23.98210.717.5−0.57.28516.6515.6323.929.92546Chhor1.910.3819.718.0125.3228.423.431.535.634.8838.244.24Chitral−4.51.67.88.56514.6225.314.224.0523.432.8338.21493Kalat−11.9−1.16.25.80812.42523.2515.1822.3522.1729.837.12007Kohat2.710.2817.917.0724.7329.314.322.330.929.5436.343.9489Mianwali0.78.617.216.8925.7330.714.524.4833.6531.593844210Multan19.619.2518.2327.0230.817.226.534.7532.5538.545.2122Muzaffarabad0.96.6513.7513.5720.7324.712.220.4829.127.5133.8339.9737Sialkot2.19.21816.6924.328.813.523.73229.7734.643.4256


### 3.1. Estimation of SPTI

At the very early stage of the computational analysis, we first prepared time series data of the SPTI index for all the stations by following the guidelines [90]. As described in Section 2.2, CDFs of the appropriate probability functions are standardized for all the stations and time scales. The minimum BIC value criterion decides the appropriateness of the probability function. Besides, we have assessed the quantile plot of theoretical and empirical densities. Thirty-two highly parameterized and extreme value distributions are included in our candidate list. The parameter estimation of these distributions and the computation of BIC values are based on Propagate [91] package of R language. Table 4 shows the BIC values calculated by fitting all the candidate distributions on SPTI-1 at selected meteorological stations, and with the lowest value of BIC, a distribution is chosen as the appropriate fitted distribution.Table 4
The candidacy of appropriate distributions according to their respective BIC values for SPTI-1.
DistributionMianwaliSialkotMuzaffarabadAstoreChitralMultanKalatChhorKohat2P beta−260.71−387.86−278.31−280.10−216.46−63.88−418.16−244.53−228.653P Weibull−472.53−718.58−555.22−480.83−400.89−317.25−638.60−444.44−622.284P beta−463.53−722.18−584.52−469.79−393.88−311.66−664.62−434.07−616.67Arcsine−211.50−322.48−364.18−272.85−194.87−38.61−353.11−178.98−216.75Burr−49.90−291.37−351.96−347.29−87.23259.98−111.45197.60−315.30Cauchy−218.97−414.95−460.90−367.96−227.51−87.98−498.98−289.68−322.84Chi−292.62−424.20−186.90−247.93−245.62113.16−220.02130.71−315.76Chi-square−321.57−546.12−456.49−441.42−271.12125.40−179.55136.71−588.65Cosine−5.68−181.37−403.03−192.31−48.26267.44−64.93206.84−159.22Curv. trapezoidal−163.25−349.64−555.47−375.44−169.25−117.58−60.78−236.04−415.05Exponential−128.46−375.68−561.96−401.81−163.06106.66−407.55−244.95−496.17F-dist−295.38−613.36−361.98−457.56−248.41158.37−180.79152.12−490.64Gamma−316.77−648.41−557.36−477.56−268.29−70.08−441.22−246.11−584.11Gen. ext. value−320.78−565.13−526.85−456.48−316.84−242.60−620.84−441.41−431.50Gen. normal−366.43−633.43−533.60−476.01−358.69−303.25−627.76−443.34−466.85Gumbel−179.72−379.06−527.12−389.62−195.10−35.02−431.30−244.43−376.46Inv. chi-square−148.95−274.23−379.11−273.05−162.6715.23−356.13−200.45−239.52Inv. Gamma−298.53−536.09−421.80−396.36−296.90−174.70−535.77−276.36−325.84Inv. Gaussian−178.04−440.80−429.50−356.23−229.22−55.60−485.95−253.33−252.24Johnson SB−364.48−633.93−547.57−478.62−357.49−298.02−623.36−445.58−463.74Johnson SU−361.60−628.62−528.86−471.37−354.24−298.52−623.36−445.58−461.94Laplace−184.73−413.96−480.10−382.24−201.82−41.51−439.14−245.85−339.58Logistic−176.53−360.75−496.70−365.07−190.99−30.47−425.80−243.25−333.03Log-normal−319.29−610.18−470.78−478.07−301.39−139.80−486.30−259.66−453.97Normal−169.86−345.46−498.92−357.27−184.07−21.42−418.91−240.70−327.01Rayleigh−171.77−356.66−524.70−373.83−186.34−23.85−420.65−241.01−358.02Scaled t-dist−241.83−423.80−494.18−369.08−242.48−122.45−535.64−356.41−330.06Skewed-normal−167.53−360.64−558.46−384.42−182.86−19.95−417.15−236.57−400.74Trapezoidal−134.44−326.85−548.58−361.71−171.10−8.78−409.37−231.50−389.97Triangular−91.67−331.66−553.32−366.35−175.4824.70−413.76−236.04−394.58Uniform−94.00−244.30−269.41−274.09−124.13170.57−80.62144.85−261.94von Mises−169.96−348.32−392.92−364.93−184.39−21.44−419.38−240.71−331.99Furthermore, in Sialkot, Muzaffarabad, and Kalat stations, the BIC values of the “four-parameter Beta” distribution are the lowest among other distributions (Sialkot, −722.18; Muzaffarabad, −584.52; Kalat, −664.62). Only in Chhor station, “Johnson SU-distribution” has given better fitness results (Chhor, −445.58). We have observed that the “Three Parameter Weibull” distribution with the lowest BIC values is the most dominant (see bold values in Table 4). For the Astore station, the histogram of the appropriate probability function, the associated quantile plot, and the temporal behavior of the standardized time series data of SPTI-1 are presented in Figure 3.Figure 3
Histogram, QQ-plot, and temporal plot of SPTI-1 at Astore station.For ease of convenience, other station plots are skipped. The red spikes indicated the drought severity and conditional dependence structure among drought episodes. The selected probability distributions with respective BIC values for all the time scales (1, 3, 6, 9, and 12) are presented in Table5. Finally, standardized time series data of SPTI for selected time scales at all the selected meteorological stations have been prepared using appropriate probability distributions.Table 5
Appropriate distributions with their respective BIC values for different time scales at nine meteorological stations.
StationsSPTI-1SPTI-3SPTI-6SPTI-9SPTI-12DistributionBICDistributionBICDistributionBICDistributionBICDistributionBICMianwali3P Weibull−472.533P Weibull−455.20Gamma−535.73Laplace−351.26Chi-square−326.56Sialkot4P beta−722.18Generalized normal−600.22Gamma−581.78Gumbel−636.77Gumbel−548.40Muzaffarabad4P beta−584.52Rayleigh−670.62Chi-square−739.36Log-normal−964.72Chi-square−866.88Astore3P Weibull−480.83Gen. extreme value−558.42Skewed-normal−745.38Rayleigh−756.25Triangular−717.03Chitral3P Weibull−400.89Generalized normal−565.65Triangular−467.11Normal−478.68Laplace−430.01Multan3P Weibull−317.254P beta−390.52Gumbel−417.03Trapezoidal−496.07Skewed-normal−355.40Kalat4P beta−664.624P beta−798.514P beta−963.704P beta−475.47Trapezoidal−500.91ChhorJohnson SU−444.443P Weibull−353.934P beta−337.40Exponential−282.64Laplace−287.91Kohat3P Weibull−622.28Gumbel−466.94Inverse Gaussian−590.98Triangular−391.63Trapezoidal−381.36


### 3.2. ELM and Its Comparative Assessment

In this research, we have assessed the performance of ELM with MLP and ARIMA in two phases. For all the individual stations selected for the current study, full-length data were divided into two independent parts, the training set and the validation set (test data). For most stations, precipitation and temperature records were available from 1951 to 2016. In the training phase, 64 years of monthly precipitation and minimum and maximum air temperature data from January 1951 to December 2013 are used to train candidate algorithms. In the testing phase, the rest of the 36 months of data (2014–2016) are considered test data to validate forecasting results. Simulations for ELM and MLP are carried out using R packagennfor [92], and forecast package for R language [93] was used to select the appropriate order of seasonal and non-seasonal specifications of the ARIMA model.Furthermore, all the ELM, MLP, and ARIMA simulations are carried out in the R 3.5.3 environment running on core i7 with a clock speed of 2.3 GHz CPU. The optimum order of ARIMA specifications with the estimated parameters is detailed in Table6. These optimum specifications were attained by running all possible ARIMA models, including all possible seasonal and non-seasonal lag values of the input time series. The ELM algorithm was trained using 23 input layer nodes and a single hidden layer with 100 hidden nodes. The algorithm is repeated 20 times, and the estimated outcomes are combined using the median operator. ELM assigns random weights to each hidden node. Furthermore, it assigns start weights to input layer nodes and generalized weights to hidden layer nodes. These weights are estimated using least absolute shrinkage and selection operator (LASSO) to keep the model parsimonious. The parametric network structure still forms a large dimension matrix, which is not feasible to illustrate numerically in a tabulated form. Furthermore, MLP is trained using two hidden layers containing 10 and 5 hidden nodes, respectively, to get optimum results.Table 6
The estimated parameters of optimum selected ARIMA models at all selected meteorological stations with (1, 3, 6, 9, and 12) month time scales.
StationsOptimum modelDriftϕ1ϕ2ϕ3ϕ4ϕ5θ1θ2θ3θ4Φ1Φ2Θ1Θ2AstoreSPTI-1ARIMA (5, 0, 1) (2, 0, 0)12—0.808−0.047−0.1230.040−0.157−0.579———0.1710.164——SPTI-3ARIMA (2, 0, 2) (2, 1, 0)12—−0.0050.114———0.8150.546——−0.611−0.315——SPTI-6ARIMA (1, 0, 0) (2, 1, 0)12−0.0020.837————————−0.648−0.350——SPTI-9ARIMA (2, 0, 1) (2, 0, 0)12—1.632−0.710———−0.576——0.2120.168——SPTI-12ARIMA (2, 0, 0) (0, 0, 1)12—1.05−0.067—————————−0.777—ChhorSPTI-1ARIMA (1, 0, 0) (2, 0, 0)12—0.179————————0.3170.219——SPTI-3ARIMA (4, 0, 0) (2, 1, 0)12—0.6960.035−0.2490.101—————−0.681−0.370——SPTI-6ARIMA (1, 0, 0) (2, 1, 0)120.0010.813————————−0.649−0.348——SPTI-9ARIMA (2, 0, 0) (2, 0, 0)12—1.000−0.191———————0.2740.162——SPTI-12ARIMA (1, 0, 1) (0, 0, 1)12—0.957————0.080—————−0.711—ChitralSPTI-1ARIMA (1, 0, 0) (1, 1, 0)12—0.123————————−0.539———SPTI-3ARIMA (2, 0, 2) (1, 1, 0)12−0.0010.0050.077———0.7820.520——−0.527———SPTI-6ARIMA (1, 0, 0) (1, 1, 0)12−0.0020.823————————−0.524———SPTI-9ARIMA (2, 0, 2) (1, 1, 0)12−0.0030.0080.788———0.9220.043——−0.540———SPTI-12ARIMA (1, 0, 0) (0, 0, 1)12—0.978——————————−0.643—KalatSPTI-1ARIMA (4, 1, 1) (0, 0, 21)12—0.3190.0380.028−0.067—−0.987—————0.2160.191SPTI-3ARIMA (2, 0, 2) (2, 1, 0)120.0010.0460.249———0.7790.449——−0.669−0.344——SPTI-6ARIMA (1, 0, 0) (2, 1, 0)12—0.880————————−0.710−0.360——SPTI-9ARIMA (0, 0, 1) (2, 0, 0)12——————0.085———0.2150.307——SPTI-12ARIMA (3, 1, 0) (0, 0, 2)12—0.1200.0770.037————————−0.782−0.052KohatSPTI-1ARIMA (1, 0, 0) (2, 0, 0)12—0.157————————0.1650.261——SPTI-3ARIMA (5, 0, 0) (2, 0, 0)12—0.912−0.104−0.5110.473−0.121————0.1410.224——SPTI-6ARIMA (4, 0, 1) (2, 0, 0)12—0.0740.838−0.096−0.209—0.959———0.1020.176——SPTI-9ARIMA (5, 0, 4) (2, 0, 0)12—−0.3400.288−0.2830.2030.4411.4491.1681.3740.9000.1300.264——SPTI-12ARIMA (1, 0, 2) (2, 0, 2)12—0.978————0.0910.067——−0.335−0.134−0.462−0.119MianwaliSPTI-1ARIMA (5, 1, 0) (2, 0, 0)12—−0.813−0.640−0.537−0.395−0.176————0.2570.304——SPTI-3ARIMA (2, 0, 3) (2, 1, 0)120.0010.5040.139———0.2540.156−0.410—−0.634−0.350——SPTI-6ARIMA (0, 1, 0) (2, 0, 0)12——————————0.2540.141——SPTI-9ARIMA (0, 1, 0) (2, 0, 0)12——————————0.3280.284——SPTI-12ARIMA (0, 1, 0) (0, 0, 1)12————————————−0.656—MultanSPTI-1ARIMA (1, 0, 2) (2, 0, 0)12—0.941————−0.861−0.116——0.2130.149——SPTI-3ARIMA (4, 0, 4) (2, 0, 0)12—−0.161—1.031−0.002−0.091—0.9411.7610.8520.6760.187———SPTI-6ARIMA (2, 0, 4) (2, 0, 0)12—0.023−0.120———0.9220.9560.8970.8610.1060.151——SPTI-9ARIMA (0, 1, 0) (2, 0, 0)12——————————0.1840.202——SPTI-12ARIMA (2, 1, 1) (1, 0, 1)12—0.251−0.05———−0.296———−0.130—−0.495—MuzaffarabadSPTI-1ARIMA (1, 0, 0) (2, 0, 0)12—0.094————————0.251———SPTI-3ARIMA (5, 0, 1) (2, 0, 0)12—1.629−0.749−0.4500.807−0.365−0.760———0.2480.241——SPTI-6ARIMA (3, 0, 1) (2, 0, 0)12—1.509−0.561−0.067——−0.553———0.2060.150——SPTI-9ARIMA (4, 0, 4) (2, 0, 0)12—0.561−0.076−0.5510.656—0.4270.4890.8230.0390.1880.324——SPTI-12ARIMA (1, 0, 2) (1, 0, 1)12—0.972————0.0360.066——−0.030—−0.690—SialkotSPTI-1ARIMA (2, 0, 2) (2, 0, 0)12—0.895−0.854———−0.7690.714——0.2310.251——SPTI-3ARIMA (3, 0, 2) (2, 1, 0)120.001−0.0240.026−0.042——0.6970.457——−0.619−0.320——SPTI-6ARIMA (2, 0, 2) (2, 0, 0)12—1.605−0.686———−0.7680.078——−0.680−0310——SPTI-9ARIMA (2, 0, 2) (1, 1, 0)120.0011.857−0.891————1.0590.180——−0.472———SPTI-12ARIMA (1, 0, 0) (0, 0, 2)12—0.968——————————−0.6650.077The MLP algorithm drastically increases as we increase the number of hidden layers or by increasing the number of nodes of hidden layers. Hence, the structure of MLP is finalized with 23 nodes of the input layer, two hidden layers with 10 and 5 nodes, respectively, and a single output layer. This parsimonious structure still forms a larger matrix of user-defined parameters. Due to complexity and numeric hazard, the estimated results for user-defined parameters are skipped. The performance of ELM, MLP, and ARIMA algorithms was assessed using performance assessment metrics, including RMSE, MAE, MAPE, KGE, and Willmott index of agreement.Table7 provides numerical results of these performance assessment metrics using training data sets for ELM, MLP, and ARIMA models at selected meteorological stations with (1, 3, 6, 9, and 12) month lead time scales. Results indicate that ELM performs better than MLP and ARIMA models. For the assessment of candidate models, the numerical results of all statistical metrics for SPTI-1 are illustrated with details. The best performance of ELM has been observed at Chitral, with a minimum value of RMSE (0.374). However, MLP best performed at Chhor station with a minimum value of RMSE as 0.598, and the ARIMA model overall best performed at the Chhor station with RMSE (0.632). While using MAE, the minimum values for ELM, MLP, and ARIMA are 0.272, 0.449, and 0.5, respectively. As for MAPE, these quantities are 259.59, 161.33, and 324.7, respectively. The values KGE for ELM, MLP, and ARIMA models at Astore station are 0.712, 0.517, and 0.314, respectively.Table 7
Prediction performance of ELM, MLP, and ARIMA algorithms using training data for SPTI with 1–12 month time scales at selected meteorological stations.
StationRMSEKGEWIMAEMAPE (%age)ELMMLPARIMAELMMLPARIMAELMMLPARIMAELMMLPARIMAELMMLPARIMASPTI-1Astore0.5230.8360.8410.7120.5170.3490.9090.7480.6640.3950.6600.671470.1519318.622633.715Chhor0.3870.5980.6320.7190.5330.3330.9100.7710.6630.2810.4490.500445.292335.141749.993Chitral0.3740.6800.7620.8380.6800.6710.9600.8470.8250.2720.5330.582259.592161.335324.699Kalat0.4020.6240.6810.7210.5560.3080.9070.7780.6430.3020.4920.558517.2201870.592903.828Kohat0.5830.9480.9670.6420.3830.1300.8800.6660.4900.4610.7660.776608.050856.133331.864Mianwali0.4160.7340.8230.7350.5030.3080.9180.7480.6250.3230.5880.679675.666978.739576.477Multan0.4890.7820.7950.6130.2920.3190.8720.6160.4250.3820.6370.674592.744365.6381005.901Muzaffarabad0.5370.7960.8370.5710.4250.1100.8910.7220.6000.4120.6360.655660.512348.687834.538Sialkot0.5230.7740.8220.7550.5920.4040.9230.7990.7080.4120.6270.662436.879238.422380.266SPTI-3Astore0.3050.5060.5580.9220.8720.8470.9780.9390.9250.2180.3770.425270.129115.508396.597Chhor0.2750.4310.4740.9050.8340.7810.9700.9230.9050.1900.3050.356266.512126.178351.577Chitral0.2520.4030.4610.9500.9170.9020.9850.9610.9510.1770.2960.34589.23646.108128.119Kalat0.3420.5040.5430.9020.8500.8240.9670.9270.9140.2350.3620.403202.44693.793164.865Kohat0.3950.5800.6250.8740.8110.7090.9590.7100.8800.2910.4310.473206.435125.313242.491Mianwali0.3300.5390.6320.9080.8330.7720.9720.9210.8890.2380.3940.472180.31976.808163.529Multan0.4160.6600.6840.8630.74606390.9520.8770.8430.2950.4850.522196.75092.480175.258Muzaffarabad0.3480.5220.5570.8250.8200.5220.9640.9170.8940.2630.3880.428249.060189.263490.751Sialkot0.3170.4760.5310.76307560.8010.9700.9270.9070.2290.3550.408230.092557.723263.315SPTI-6Astore0.2870.3510.4550.9450.9120.8940.9910.9810.9490.1980.2480.34775.91886.424129.332Chhor0.2550.3530.4300.9390.9160.8590.9790.9600.9360.1550.2290.30798.630131.632181.744Chitral0.2470.2990.4020.9640.9550.9250.9860.9780.9620.1670.2110.30681.65152.158131.476Kalat0.3160.3800.4390.9560.9420.8970.9930.9760.9520.1930.2400.30179.513300.351238.315Kohat0.8830.4550.5070.9220.9570.8230.9680.9830.9300.6950.3130.370173.514253.266703.016Mianwali0.3190.4090.4880.9120.8980.8660.9720.9530.9330.2180.2720.332209.50498.576422.293Multan0.3820.4780.4810.9220.8720.8410.9930.9630.9370.2640.3180.337179.67699.793136.787Muzaffarabad0.3900.4710.5270.8850.7530.3430.9590.9440.9140.2780.3370.395158.979114.439364.869Sialkot0.2980.3750.4320.8890.8820.6910.9750.9590.9440.1990.2490.314117.07897.213149.507SPTI-9Astore0.2880.3460.4340.9490.9360.8740.9800.9710.9510.1970.2480.30976.93362.029194.755Chhor0.2830.3550.4330.9430.9100.7640.9770.9510.9070.1580.2210.289150.82892.464129.091Chitral0.2630.3080.4430.9580.9550.9140.9840.9790.9570.1850.2220.319116.52758.64586.048Kalat0.2840.3440.4090.9510.9450.9270.9820.9740.9630.1690.2100.241115.838171.454201.090Kohat0.3180.3710.4150.9420.7650.8460.9750.8750.9390.2200.2590.30289.845123.001117.730Mianwali0.2640.3070.3880.9480.9480.9160.9820.9750.9610.1840.2190.277466.11896.627164.084Multan0.3660.4360.5470.9280.9050.8620.9680.9530.9280.2470.2940.352143.63490.969124.937Muzaffarabad0.3230.3850.4710.9380.9210.8490.9750.9640.9410.2400.2890.357180.985107.083712.336Sialkot0.2980.3340.4230.9470.9380.8950.9760.9700.9510.2020.2340.310121.5391223.05552.815SPTI-12Astore0.3110.3360.3000.9510.9450.9380.9730.9770.9780.2030.2160.204115.137158.271131.494Chhor0.3020.3400.3020.9460.9330.8640.9740.9670.9730.1510.1540.15897.432241.04060.678Chitral0.3480.3720.3400.9400.9360.9230.9710.9680.9730.2140.2180.21394.94283.428144.242Kalat0.2400.2840.2570.9730.9680.9660.9890.9850.9870.1220.1310.12451.69059.213457.499Kohat0.3010.3420.3020.9520.9440.9390.9780.9720.9780.2030.2240.204146.08178.430131.950Mianwali0.2540.2770.2540.9680.9590.3280.9850.9820.9840.1690.1800.16797.51894.95070.970Multan0.3420.3830.3620.9350.8140.8990.9570.8130.9620.2160.2290.225412.64480.412142.043Muzaffarabad0.3440.3700.3350.9410.9350.9270.9720.9680.9720.2460.2600.242111.26090.649165.662Sialkot0.3290.3640.3290.9380.9320.9110.9720.9660.9710.2090.2240.213108.491114.55492.435Similarly, the numeric quantities of WI at the Astore station are 0.999, 0.748, and 0.664, respectively. The KGE index indorses ELM’s superior performance at all the stations by providing maximum values as compared to MLP and ARIMA. Similarly, Willmott’s agreement “WI” index consistently provides the highest values for the ELM algorithm. KGE and WI are considered the most appropriate metrics for the performance assessment of meteorological and hydrological models. The KGE is calculated using the correlation coefficient, the ratio of variations, and the ratio of averages of predicted and observed series using equation (9). The values of “WI” for the ELM model for all the selected stations are close to 1, which endorses the ELM as the best performing model. The similar superior performance of ELM continued for other time scales at selected stations. Overall results for the training phase show that the ELM model has shown good agreement at all selected stations.As the time scale of the drought index increases, the performance of the proposed algorithm improves. As a result, the ELM algorithm showed superior performance to its competitive algorithms (MLP and ARIMA). A comparison of all the performance assessment metrics concluded that ELM algorithm is selected as the adequate model for the estimation and forecasting of drought indices (see Table7).The consistency and co-movement of the observed and estimated values of SPTI are further assessed while employing Karl Pearson’s product-moment correlation coefficient. Table8 shows the numerical results of the correlation coefficient between the observed and predicted values of SPTI using ELM, MLP, and ARIMA models for training data. The quantitative results of the Astore station for ELM, MLP, and ARIMA are 0.87, 0.59, and 0.56, respectively, indicating a better agreement of ELM to predict SPTI-1 contrary to other candidate algorithms. At Chitral station, values of correlation are (0.93, 0.74, and 0.69). For SPTI-3 at Astore station, these results are 0.96, 0.88, and 0.86, respectively, and at Chitral station, the correlation values for ELM, MLP, and ARIMA are 0.97, 0.93, and 0.90. These numerical results clearly illustrated that the performance of the ELM model significantly improved as the time scale increased. A similar pattern of the superior prediction performance of ELM continued for other time scales at all the selected stations. Although all the models have shown reasonable prediction performance, quantitative results evidently confirmed that the estimated values of SPTI using the ELM model strongly correlate with the observed values of SPTI for any time scale.Table 8
Pearson’s correlation coefficient between the observed and predicted values of SPTI 1–12 month time scales using ELM, MLP, and ARIMA algorithms at selected meteorological stations.
ScaleAstoreChhorChitralKalatKohatMianwaliMultanMuzaffarabadSialkotELMMLPARIMAELMMLPARIMAELMMLPARIMAELMMLPARIMAELMMLPARIMAELMMLPARIMAELMMLPARIMAELMMLPARIMAELMMLPARIMASPTI-10.870.590.560.860.620.540.930.740.690.860.630.520.840.480.390.880.590.430.840.410.320.840.560.480.880.660.60SPTI-30.960.880.860.950.860.830.970.930.900.940.860.840.920.830.800.950.850.800.910.780.750.930.850.820.940.860.83SPTI-60.960.940.900.960.920.880.970.960.930.950.930.910.940.900.870.940.910.870.930.890.890.930.890.860.950.920.89SPTI-90.960.940.910.960.280.840.970.960.920.970.950.930.950.770.890.960.950.920.940.910.860.950.930.890.950.940.91SPTI-120.960.950.960.950.930.950.940.940.950.980.970.980.960.950.960.970.960.970.770.920.930.950.940.950.950.930.94Usually, climatic and meteorological studies comprise high-frequency datasets that require fast algorithms. So, speed is a notable characteristic for determining the reliability of the algorithm. The algorithm selection for climatic studies is subjective in terms of speed and relative efficiency. ELM has the novelty of being the fastest algorithm among the ANN class to solve complex datasets. ELM algorithm training and testing time were almost 32 times faster than ANN, indicating ELM's supremacy over other ANN algorithms [58].The functional relationship between the actual (observed) and predicted values of SPTI using ELM and other algorithms for the “Astore” station is shown in Figure4 using a line graph, which depicts significantly less variation among the observed values of SPTI and the predicted values using ELM.Figure 4
Performance assessment of ELM, MLP, and ARIMA models through co-movement of observed and predicted values of SPTI at Astore station for the training phase.MLP and ARIMA were unable to capture all the shocks in historical values of SPTI, and departure from observed values was significant. Here, the ELM algorithm reflects more precise and accurate predictions. Figures5– detail the ELM's prediction performance for all the time scales at the selected stations, which depicts the ELM algorithm's better prediction performance.Figure 5
Co-movement of the observed and predicted values of SPTI using ELM at Astore, Chhor, and Chitral stations for the training phase.Figure 6
Co-movement of the observed and predicted values of SPTI using ELM at Kalat, Kohat, and Mianwali stations for the training phase.Figure 7
Co-movement of the observed and predicted values of SPTI using ELM at Muzaffarabad, Multan, and Sialkot stations for the training phase.It was observed that as the time scale increases, the prediction performance of algorithms substantially improves. For SPTI-12, significantly fewer deviations have been observed among drought index’s paired (observed and predicted) values. These multi-line plots indicate that the ELM model incurred smaller errors than the two counterparts. The predicted values using ELM model follow the observed values of SPTI more precisely. These graphical presentations show that the ELM model has shown more accuracy than MLP and ARIMA models for estimating SPTI for various lead time scales at selected meteorological stations.Scatter plots of the observed and predicted time series data sets are another way to assess the prediction performance of probabilistic, machine learning, and ANN algorithms. Figures8–10 show the scatter plots of the historical observed and predicted values of SPTI-6 using ELM, MLP, and ARIMA models at all the selected meteorological stations. Predicted values through the ELM model have shown a strong correlation with the observed values of SPTI through scatter plots. Another graphical presentation endorses the superior performance of ELM. We can observe that the ELM algorithm showed more accuracy and can potentially predict drought conditions in any climatic zone.Figure 8
Scatter plots of the observed and predicted values of SPTI-6 using ELM, MLP, and ARIMA at Astore, Chhor, and Chitral stations for the training phase.Figure 9
Scatter plots of the observed and predicted values of SPTI using ELM, MLP, and ARIMA at Kalat, Kohat, and Mianwali stations for the training phase.Figure 10
Scatter plots of the observed and predicted values of SPTI using ELM, MLP, and ARIMA at Muzaffarabad, Multan, and Sialkot stations for the training phase.Figure11 represents the Taylor diagrams for the Astore station with (1, 3, 6, 9, and 12) month time scales for the training phase. Taylor diagrams are the more comprehensive and precise way to represent the estimation and forecast ability of a model. The similarity between the predicted and actual values of SPTI is evaluated in terms of correlation (as a measure of closeness), and the variation is assessed by the standard deviation (SD) and the RMSE. For SPTI-1, the correlation of modeled data using the ELM algorithm with actual observations was about 0.9, followed by MLP and ARIMA with 0.6 each. As time scale increases, the prediction performance of algorithms significantly improves. For SPTI-3, the ELM algorithm is significantly closer to the actual values as its correlation is about (0.97) as compared to MLP (0.9) and ARIMA (0.85). The Taylor diagram exhibits the superior performance of ELM algorithm for estimating SPTI with (1, 3, 6, 9, and 12) month lead time scales (see Figure 11). Figure 12 illustrates the violin plots related to the training phase of ELM and other candidate models. The red dot represents the mean value, a thick white bar represents the interquartile, and a thin blue line represents the whole data set distribution. These are the components of the boxplot, but the colored area on both sides of the blue line is the estimated kernel density of the shape of the distribution of data. Higher probabilities are associated with wider parts of the diagram, and thinner sections show lower probabilities. This is another systematic way to compare the prediction performance of models. It is noticeable that the mean of the predicted and observed values of SPTI-6 for the ELM model was almost similar for all the stations, whereas other models have slight variations. All the models have shown reasonable agreement but Violin diagrams of ELM and observed data are nearly identical. These graphical illustrations affirmed that the proposed ELM model is better in estimating the actual values of SPTI at all selected meteorological stations.Figure 11
Taylor diagrams for the evaluation of the prediction performance of ELM, MLP, and ARIMA algorithms at Astore station with different time scales in the training phase.Figure 12
Violin plots of the observed SPTI-6 and the modeled data using ELM, MLP, and ARIMA models at all the selected meteorological stations in the training phase.After calibrating and validating algorithms for training datasets, the proposed algorithm's generalization capability has been assessed in the next step. Finally, out of the sample forecast of SPTI for all the lead time scales, (1–12) is carried out for 36 months from 2014 to 2016. These forecasts are considered sufficient for drought preparedness and mitigation policies. Similar performance metrics have been used to analyze the difference between the observed and forecasted values of SPTI for different time scales at all selected meteorological stations. Numerical results related to these performance metrics are given in Table9. If the comparison of KGE was made at Astore station, the values of KGE for ELM, MLP, and ARIMA models are 0.988, 0.771, and 0.671, respectively. The KGE for the ELM model is significantly higher than its counterparts for all the time scales at all the selected meteorological stations, which clearly endorses the better forecast capability of the ELM model. Similarly, WI and MAE quantitative results also ratified that the ELM model outperformed the ANN and ARIMA models. RMSE and MAPE endorsed ELM as a better forecasting model for most stations.Table 9
Forecast performance of ELM, MLP, and ARIMA algorithms for SPTI with 1–12 month time scales at nine meteorological stations.
StationRMSEKGEWIMAEMAPE (%age)ELMMLPARIMAELMMLPARIMAELMMLPARIMAELMMLPARIMAELMMLPARIMASPTI-1Astore0.6881.0340.8930.9880.7710.6710.9970.8970.8840.7980.8550.865247.353251.284364.431Chhor0.5790.6000.5590.8710.7120.5410.9990.7300.6630.5310.4470.468559.287269.288396.060Chitral0.6430.6470.6900.9090.6850.6780.9790.8360.8230.4950.5100.520142.102193.952414.532Kalat0.7760.8320.7270.9520.7790.7670.9980.9660.9760.6760.7310.684500.369824.503616.762Kohat0.9560.9701.0310.9790.9300.9000.9990.9920.9540.7910.8150.869528.0661136.164268.277Mianwali0.7730.9840.8010.9610.7530.6900.9980.9540.9230.6520.8180.655309.146369.444860.231Multan0.7140.6660.8240.9860.4280.6870.9970.7400.8850.6120.5770.706545.5681291.2043958.975Muzaffarabad0.7640.7430.7240.9370.8770.6730.9970.9660.9940.6050.5970.553321.8103469.574455.340Sialkot0.7200.7450.8220.9310.9020.9140.9950.9890.9820.5540.5930.66391.037151.493269.868SPTI-3Astore0.7640.8960.6070.9700.7440.5110.9910.8650.7600.5890.7350.488797.80063.99284.646Chhor0.5630.5030.5390.86609620.6780.8210.8730.8710.4380.3660.385193.072159.359141.894Chitral0.4680.5220.5140.9210.7090.8200.9650.9080.9180.4080.3960.4052354.0679.394150.935Kalat0.9871.2471.1520.9960.8580.6780.9990.9560.7580.8501.0670.899329.382314.217199.070Kohat0.7710.7750.8730.9740.9650.730.9940.9930.8910.7290.6680.726116.387330.497134.998Mianwali1.1041.2211.0410.9970.9610.9430.9970.9870.9560.9141.0140.8061082.42561.266173.703Multan0.8010.7900.9820.8940.9270.9380.9990.9970.9780.6600.6540.826210.083193.947210.070Muzaffarabad0.7910.7850.7770.9670.8560.8730.9990.9900.9960.6020.6860.601230.935809.403334.650Sialkot0.6810.7250.7680.9380.9240.6280.9960.9770.8110.5200.5320.601265.276259.027201.413SPTI-6Astore0.4020.6700.4820.9620.6070.5410.9850.7820.7330.4920.5820.40846.301343.20737.464Chhor0.8190.8890.6690.9100.6870.4560.9980.9670.8760.7020.7710.475697.320269.54173.976Chitral0.3740.4050.3990.9610.8620.8460.9840.9350.9430.3080.3250.308172.828140.348140.661Kalat0.9720.7921.2310.9650.9210.9330.9930.8840.9760.8070.6661.0221364.65335.608184.772Kohat0.9570.9300.8690.9240.8760.7590.9900.9480.9180.7640.8410.71989.183301.444109.716Mianwali1.4112.1481.4570.9740.7880.8850.9990.9750.9871.3191.8681.145347.413215.620298.005Multan0.7831.2091.0760.9310.7810.9260.9990.9990.9910.6730.9820.886222.776194.657173.698Muzaffarabad0.7711.1760.8720.9890.9350.8930.9990.9950.9810.6480.8220.6952458.04214.1594533.34Sialkot0.8800.7330.6560.9530.9440.9280.9980.9910.9940.7560.6700.557776.441532.295178.914SPTI-9Astore0.6340.6470.8740.9680.7280.7830.9900.8570.8690.4960.5420.77641.23049.977101.519Chhor1.1341.3280.9610.9490.9360.9450.9940.9960.8290.9441.1510.760358.134243.020584.953Chitral0.3720.5040.3830.9530.8630.8200.9840.9210.9100.3220.4060.29753.6963157.8534.444Kalat0.4281.1330.7820.9910.7240.8450.9990.8920.9350.8760.9820.713187.052360.12550.731Kohat1.0321.0420.8360.9410.9320.7740.9920.9870.9620.8400.7910.61798.19087.81372.893Mianwali1.3502.0621.4570.9700.9240.6160.9980.9820.9641.4291.8951.145253.908162.066298.005Multan0.7460.6901.1480.9730.8730.8980.9990.9980.9990.6090.5950.9461139.221061.42134.464Muzaffarabad0.7630.8000.7610.9850.9480.9380.9990.9790.9900.6340.6600.59278.798144.12688.432Sialkot0.7650.6810.8390.9250.9040.8510.9970.9960.9850.5890.4850.65594.806227.032221.098SPTI-12Astore0.2580.7541.2290.9620.7210.7700.9920.9060.8690.2150.6581.20210.93668.122187.249Chhor0.6700.7070.7280.9870.9750.8790.9960.9350.9210.5660.6000.648197.9181307.502828.820Chitral0.5210.5240.7720.9960.7470.6780.9990.9220.8540.4060.4130.65338.71035.631117.257Kalat1.3461.2161.1680.9750.9940.6120.9980.9980.8041.1791.0260.9841739.84558.414343.401Kohat0.6151.0660.9630.9810.9590.9000.9970.9950.9810.8600.7900.81984.52266.743144.398Mianwali1.4021.5131.4840.9940.8820.7750.9990.9890.8371.6291.2981.273147.015180.588169.143Multan0.5611.0720.6640.9950.9160.8700.9980.9740.9700.4340.9970.55360.9761233.10223.381Muzaffarabad0.7201.1170.7540.9850.7560.8900.9970.8990.9620.5640.7800.58649.482257.15694.969Sialkot0.7930.7270.7660.9450.9260.9220.9930.9920.9730.6300.5930.59361.040181.0031734.785The functional relationship between the observed and forecasted values of SPTI using ELM, MLP, and ARIMA models for all the time scales at Astore stations is illustrated in Figure13 for the test phase starting from January 2014 to December 2016 for 36 months. A significantly smaller degree of deviation in SPTI for the ELM model and observed values of SPTI were exhibited. Although MLP has shown a reasonably good forecast performance compared to the stochastic seasonal ARIMA model, the ELM model evidently shows superior forecast performance.Figure 13
Observed and forecasted SPTI by ELM, MLP, and ARIMA models at Astore station for 1, 3, 6, 9, and 12 month lead time scales.In order to check the appropriateness of the ELM model, scatter plots were prepared using time series data of the observed and forecasted SPTI at Astore station for all lead time scales (Figure14). These scatter plots have shown a significant difference in the forecast performance of ELM and other models. The scatter plots depict the correlation, goodness-of-fit, and the extent of agreement between the observed and forecasted SPTI. The ELM model also clearly outperformed MLP and ARIMA models for the testing phase for all the selected time scales.Figure 14
Scatter plots of the observed and forecasted values of SPTI using ELM, MLP, and ARIMA models at Astore station with all lead time scales for the test phase. The least-squares regression line is also included.


### 3.3. Discussion

Drought is a multifaceted and commonly occurring hazard in several parts of the world. Its impacts are prevalent in the agriculture, socio-economic, and energy sectors. However, precise drought monitoring and estimation techniques can assist in decreasing the vulnerability of society to drought. The primary objective of the current study was to test the appropriateness and usefulness of the ELM model relative to other ANN (MLP) models and stochastic (ARIMA) models for drought forecasting. The prediction and forecast performance of ELM is compared with other ANN algorithms (MLP) and statistical stochastic (ARIMA) models. The prediction and forecast performance of models is assessed using numerous performance metrics, including RMSE, MAE, MAPE, KGE, WI, and Karl Pearson's correlation coefficient. The quantitative assessment revealed that both the ANN models (ELM and MLP) performed better than the stochastic model (ARIMA), and among the ANN models, ELM has shown supremacy by producing the smallest RMSE, MAE, and MAPE values and the maximum values for KGE, WI, and correlation coefficient for almost all the meteorological stations.Furthermore, ELM shows better agreement for both the training and test phases to predict the SPTI at all climatic stations than its counterparts. The efficiency of ELM, contrary to other models, is evident based on the performance metrics. A similar forecast performance has continued for higher-order time scales, consistent with earlier studies [30, 58]. Computational time consumed by drought modeling algorithms also needs to be optimized. Usually, large datasets are used as input variables for real-time drought modeling, which affects the computational performance of different models in terms of time. The ELM model is significantly faster than its counterparts. By evaluating all the numerical results of performance metrics and different graphical illustrations, it can be easily concluded that the ELM model attains the most accurate drought forecasting performance during training and test phases. The study revealed that ELM is the most appropriate, reliable, and efficient algorithm for drought prediction and forecasting. This study suggests that ELM can be used as an early warning drought forecasting tool for developing drought mitigation policies.


## 3.1. Estimation of SPTI

At the very early stage of the computational analysis, we first prepared time series data of the SPTI index for all the stations by following the guidelines [90]. As described in Section 2.2, CDFs of the appropriate probability functions are standardized for all the stations and time scales. The minimum BIC value criterion decides the appropriateness of the probability function. Besides, we have assessed the quantile plot of theoretical and empirical densities. Thirty-two highly parameterized and extreme value distributions are included in our candidate list. The parameter estimation of these distributions and the computation of BIC values are based on Propagate [91] package of R language. Table 4 shows the BIC values calculated by fitting all the candidate distributions on SPTI-1 at selected meteorological stations, and with the lowest value of BIC, a distribution is chosen as the appropriate fitted distribution.Table 4
The candidacy of appropriate distributions according to their respective BIC values for SPTI-1.
DistributionMianwaliSialkotMuzaffarabadAstoreChitralMultanKalatChhorKohat2P beta−260.71−387.86−278.31−280.10−216.46−63.88−418.16−244.53−228.653P Weibull−472.53−718.58−555.22−480.83−400.89−317.25−638.60−444.44−622.284P beta−463.53−722.18−584.52−469.79−393.88−311.66−664.62−434.07−616.67Arcsine−211.50−322.48−364.18−272.85−194.87−38.61−353.11−178.98−216.75Burr−49.90−291.37−351.96−347.29−87.23259.98−111.45197.60−315.30Cauchy−218.97−414.95−460.90−367.96−227.51−87.98−498.98−289.68−322.84Chi−292.62−424.20−186.90−247.93−245.62113.16−220.02130.71−315.76Chi-square−321.57−546.12−456.49−441.42−271.12125.40−179.55136.71−588.65Cosine−5.68−181.37−403.03−192.31−48.26267.44−64.93206.84−159.22Curv. trapezoidal−163.25−349.64−555.47−375.44−169.25−117.58−60.78−236.04−415.05Exponential−128.46−375.68−561.96−401.81−163.06106.66−407.55−244.95−496.17F-dist−295.38−613.36−361.98−457.56−248.41158.37−180.79152.12−490.64Gamma−316.77−648.41−557.36−477.56−268.29−70.08−441.22−246.11−584.11Gen. ext. value−320.78−565.13−526.85−456.48−316.84−242.60−620.84−441.41−431.50Gen. normal−366.43−633.43−533.60−476.01−358.69−303.25−627.76−443.34−466.85Gumbel−179.72−379.06−527.12−389.62−195.10−35.02−431.30−244.43−376.46Inv. chi-square−148.95−274.23−379.11−273.05−162.6715.23−356.13−200.45−239.52Inv. Gamma−298.53−536.09−421.80−396.36−296.90−174.70−535.77−276.36−325.84Inv. Gaussian−178.04−440.80−429.50−356.23−229.22−55.60−485.95−253.33−252.24Johnson SB−364.48−633.93−547.57−478.62−357.49−298.02−623.36−445.58−463.74Johnson SU−361.60−628.62−528.86−471.37−354.24−298.52−623.36−445.58−461.94Laplace−184.73−413.96−480.10−382.24−201.82−41.51−439.14−245.85−339.58Logistic−176.53−360.75−496.70−365.07−190.99−30.47−425.80−243.25−333.03Log-normal−319.29−610.18−470.78−478.07−301.39−139.80−486.30−259.66−453.97Normal−169.86−345.46−498.92−357.27−184.07−21.42−418.91−240.70−327.01Rayleigh−171.77−356.66−524.70−373.83−186.34−23.85−420.65−241.01−358.02Scaled t-dist−241.83−423.80−494.18−369.08−242.48−122.45−535.64−356.41−330.06Skewed-normal−167.53−360.64−558.46−384.42−182.86−19.95−417.15−236.57−400.74Trapezoidal−134.44−326.85−548.58−361.71−171.10−8.78−409.37−231.50−389.97Triangular−91.67−331.66−553.32−366.35−175.4824.70−413.76−236.04−394.58Uniform−94.00−244.30−269.41−274.09−124.13170.57−80.62144.85−261.94von Mises−169.96−348.32−392.92−364.93−184.39−21.44−419.38−240.71−331.99Furthermore, in Sialkot, Muzaffarabad, and Kalat stations, the BIC values of the “four-parameter Beta” distribution are the lowest among other distributions (Sialkot, −722.18; Muzaffarabad, −584.52; Kalat, −664.62). Only in Chhor station, “Johnson SU-distribution” has given better fitness results (Chhor, −445.58). We have observed that the “Three Parameter Weibull” distribution with the lowest BIC values is the most dominant (see bold values in Table 4). For the Astore station, the histogram of the appropriate probability function, the associated quantile plot, and the temporal behavior of the standardized time series data of SPTI-1 are presented in Figure 3.Figure 3
Histogram, QQ-plot, and temporal plot of SPTI-1 at Astore station.For ease of convenience, other station plots are skipped. The red spikes indicated the drought severity and conditional dependence structure among drought episodes. The selected probability distributions with respective BIC values for all the time scales (1, 3, 6, 9, and 12) are presented in Table5. Finally, standardized time series data of SPTI for selected time scales at all the selected meteorological stations have been prepared using appropriate probability distributions.Table 5
Appropriate distributions with their respective BIC values for different time scales at nine meteorological stations.
StationsSPTI-1SPTI-3SPTI-6SPTI-9SPTI-12DistributionBICDistributionBICDistributionBICDistributionBICDistributionBICMianwali3P Weibull−472.533P Weibull−455.20Gamma−535.73Laplace−351.26Chi-square−326.56Sialkot4P beta−722.18Generalized normal−600.22Gamma−581.78Gumbel−636.77Gumbel−548.40Muzaffarabad4P beta−584.52Rayleigh−670.62Chi-square−739.36Log-normal−964.72Chi-square−866.88Astore3P Weibull−480.83Gen. extreme value−558.42Skewed-normal−745.38Rayleigh−756.25Triangular−717.03Chitral3P Weibull−400.89Generalized normal−565.65Triangular−467.11Normal−478.68Laplace−430.01Multan3P Weibull−317.254P beta−390.52Gumbel−417.03Trapezoidal−496.07Skewed-normal−355.40Kalat4P beta−664.624P beta−798.514P beta−963.704P beta−475.47Trapezoidal−500.91ChhorJohnson SU−444.443P Weibull−353.934P beta−337.40Exponential−282.64Laplace−287.91Kohat3P Weibull−622.28Gumbel−466.94Inverse Gaussian−590.98Triangular−391.63Trapezoidal−381.36


## 3.2. ELM and Its Comparative Assessment

In this research, we have assessed the performance of ELM with MLP and ARIMA in two phases. For all the individual stations selected for the current study, full-length data were divided into two independent parts, the training set and the validation set (test data). For most stations, precipitation and temperature records were available from 1951 to 2016. In the training phase, 64 years of monthly precipitation and minimum and maximum air temperature data from January 1951 to December 2013 are used to train candidate algorithms. In the testing phase, the rest of the 36 months of data (2014–2016) are considered test data to validate forecasting results. Simulations for ELM and MLP are carried out using R packagennfor [92], and forecast package for R language [93] was used to select the appropriate order of seasonal and non-seasonal specifications of the ARIMA model.Furthermore, all the ELM, MLP, and ARIMA simulations are carried out in the R 3.5.3 environment running on core i7 with a clock speed of 2.3 GHz CPU. The optimum order of ARIMA specifications with the estimated parameters is detailed in Table6. These optimum specifications were attained by running all possible ARIMA models, including all possible seasonal and non-seasonal lag values of the input time series. The ELM algorithm was trained using 23 input layer nodes and a single hidden layer with 100 hidden nodes. The algorithm is repeated 20 times, and the estimated outcomes are combined using the median operator. ELM assigns random weights to each hidden node. Furthermore, it assigns start weights to input layer nodes and generalized weights to hidden layer nodes. These weights are estimated using least absolute shrinkage and selection operator (LASSO) to keep the model parsimonious. The parametric network structure still forms a large dimension matrix, which is not feasible to illustrate numerically in a tabulated form. Furthermore, MLP is trained using two hidden layers containing 10 and 5 hidden nodes, respectively, to get optimum results.Table 6
The estimated parameters of optimum selected ARIMA models at all selected meteorological stations with (1, 3, 6, 9, and 12) month time scales.
StationsOptimum modelDriftϕ1ϕ2ϕ3ϕ4ϕ5θ1θ2θ3θ4Φ1Φ2Θ1Θ2AstoreSPTI-1ARIMA (5, 0, 1) (2, 0, 0)12—0.808−0.047−0.1230.040−0.157−0.579———0.1710.164——SPTI-3ARIMA (2, 0, 2) (2, 1, 0)12—−0.0050.114———0.8150.546——−0.611−0.315——SPTI-6ARIMA (1, 0, 0) (2, 1, 0)12−0.0020.837————————−0.648−0.350——SPTI-9ARIMA (2, 0, 1) (2, 0, 0)12—1.632−0.710———−0.576——0.2120.168——SPTI-12ARIMA (2, 0, 0) (0, 0, 1)12—1.05−0.067—————————−0.777—ChhorSPTI-1ARIMA (1, 0, 0) (2, 0, 0)12—0.179————————0.3170.219——SPTI-3ARIMA (4, 0, 0) (2, 1, 0)12—0.6960.035−0.2490.101—————−0.681−0.370——SPTI-6ARIMA (1, 0, 0) (2, 1, 0)120.0010.813————————−0.649−0.348——SPTI-9ARIMA (2, 0, 0) (2, 0, 0)12—1.000−0.191———————0.2740.162——SPTI-12ARIMA (1, 0, 1) (0, 0, 1)12—0.957————0.080—————−0.711—ChitralSPTI-1ARIMA (1, 0, 0) (1, 1, 0)12—0.123————————−0.539———SPTI-3ARIMA (2, 0, 2) (1, 1, 0)12−0.0010.0050.077———0.7820.520——−0.527———SPTI-6ARIMA (1, 0, 0) (1, 1, 0)12−0.0020.823————————−0.524———SPTI-9ARIMA (2, 0, 2) (1, 1, 0)12−0.0030.0080.788———0.9220.043——−0.540———SPTI-12ARIMA (1, 0, 0) (0, 0, 1)12—0.978——————————−0.643—KalatSPTI-1ARIMA (4, 1, 1) (0, 0, 21)12—0.3190.0380.028−0.067—−0.987—————0.2160.191SPTI-3ARIMA (2, 0, 2) (2, 1, 0)120.0010.0460.249———0.7790.449——−0.669−0.344——SPTI-6ARIMA (1, 0, 0) (2, 1, 0)12—0.880————————−0.710−0.360——SPTI-9ARIMA (0, 0, 1) (2, 0, 0)12——————0.085———0.2150.307——SPTI-12ARIMA (3, 1, 0) (0, 0, 2)12—0.1200.0770.037————————−0.782−0.052KohatSPTI-1ARIMA (1, 0, 0) (2, 0, 0)12—0.157————————0.1650.261——SPTI-3ARIMA (5, 0, 0) (2, 0, 0)12—0.912−0.104−0.5110.473−0.121————0.1410.224——SPTI-6ARIMA (4, 0, 1) (2, 0, 0)12—0.0740.838−0.096−0.209—0.959———0.1020.176——SPTI-9ARIMA (5, 0, 4) (2, 0, 0)12—−0.3400.288−0.2830.2030.4411.4491.1681.3740.9000.1300.264——SPTI-12ARIMA (1, 0, 2) (2, 0, 2)12—0.978————0.0910.067——−0.335−0.134−0.462−0.119MianwaliSPTI-1ARIMA (5, 1, 0) (2, 0, 0)12—−0.813−0.640−0.537−0.395−0.176————0.2570.304——SPTI-3ARIMA (2, 0, 3) (2, 1, 0)120.0010.5040.139———0.2540.156−0.410—−0.634−0.350——SPTI-6ARIMA (0, 1, 0) (2, 0, 0)12——————————0.2540.141——SPTI-9ARIMA (0, 1, 0) (2, 0, 0)12——————————0.3280.284——SPTI-12ARIMA (0, 1, 0) (0, 0, 1)12————————————−0.656—MultanSPTI-1ARIMA (1, 0, 2) (2, 0, 0)12—0.941————−0.861−0.116——0.2130.149——SPTI-3ARIMA (4, 0, 4) (2, 0, 0)12—−0.161—1.031−0.002−0.091—0.9411.7610.8520.6760.187———SPTI-6ARIMA (2, 0, 4) (2, 0, 0)12—0.023−0.120———0.9220.9560.8970.8610.1060.151——SPTI-9ARIMA (0, 1, 0) (2, 0, 0)12——————————0.1840.202——SPTI-12ARIMA (2, 1, 1) (1, 0, 1)12—0.251−0.05———−0.296———−0.130—−0.495—MuzaffarabadSPTI-1ARIMA (1, 0, 0) (2, 0, 0)12—0.094————————0.251———SPTI-3ARIMA (5, 0, 1) (2, 0, 0)12—1.629−0.749−0.4500.807−0.365−0.760———0.2480.241——SPTI-6ARIMA (3, 0, 1) (2, 0, 0)12—1.509−0.561−0.067——−0.553———0.2060.150——SPTI-9ARIMA (4, 0, 4) (2, 0, 0)12—0.561−0.076−0.5510.656—0.4270.4890.8230.0390.1880.324——SPTI-12ARIMA (1, 0, 2) (1, 0, 1)12—0.972————0.0360.066——−0.030—−0.690—SialkotSPTI-1ARIMA (2, 0, 2) (2, 0, 0)12—0.895−0.854———−0.7690.714——0.2310.251——SPTI-3ARIMA (3, 0, 2) (2, 1, 0)120.001−0.0240.026−0.042——0.6970.457——−0.619−0.320——SPTI-6ARIMA (2, 0, 2) (2, 0, 0)12—1.605−0.686———−0.7680.078——−0.680−0310——SPTI-9ARIMA (2, 0, 2) (1, 1, 0)120.0011.857−0.891————1.0590.180——−0.472———SPTI-12ARIMA (1, 0, 0) (0, 0, 2)12—0.968——————————−0.6650.077The MLP algorithm drastically increases as we increase the number of hidden layers or by increasing the number of nodes of hidden layers. Hence, the structure of MLP is finalized with 23 nodes of the input layer, two hidden layers with 10 and 5 nodes, respectively, and a single output layer. This parsimonious structure still forms a larger matrix of user-defined parameters. Due to complexity and numeric hazard, the estimated results for user-defined parameters are skipped. The performance of ELM, MLP, and ARIMA algorithms was assessed using performance assessment metrics, including RMSE, MAE, MAPE, KGE, and Willmott index of agreement.Table7 provides numerical results of these performance assessment metrics using training data sets for ELM, MLP, and ARIMA models at selected meteorological stations with (1, 3, 6, 9, and 12) month lead time scales. Results indicate that ELM performs better than MLP and ARIMA models. For the assessment of candidate models, the numerical results of all statistical metrics for SPTI-1 are illustrated with details. The best performance of ELM has been observed at Chitral, with a minimum value of RMSE (0.374). However, MLP best performed at Chhor station with a minimum value of RMSE as 0.598, and the ARIMA model overall best performed at the Chhor station with RMSE (0.632). While using MAE, the minimum values for ELM, MLP, and ARIMA are 0.272, 0.449, and 0.5, respectively. As for MAPE, these quantities are 259.59, 161.33, and 324.7, respectively. The values KGE for ELM, MLP, and ARIMA models at Astore station are 0.712, 0.517, and 0.314, respectively.Table 7
Prediction performance of ELM, MLP, and ARIMA algorithms using training data for SPTI with 1–12 month time scales at selected meteorological stations.
StationRMSEKGEWIMAEMAPE (%age)ELMMLPARIMAELMMLPARIMAELMMLPARIMAELMMLPARIMAELMMLPARIMASPTI-1Astore0.5230.8360.8410.7120.5170.3490.9090.7480.6640.3950.6600.671470.1519318.622633.715Chhor0.3870.5980.6320.7190.5330.3330.9100.7710.6630.2810.4490.500445.292335.141749.993Chitral0.3740.6800.7620.8380.6800.6710.9600.8470.8250.2720.5330.582259.592161.335324.699Kalat0.4020.6240.6810.7210.5560.3080.9070.7780.6430.3020.4920.558517.2201870.592903.828Kohat0.5830.9480.9670.6420.3830.1300.8800.6660.4900.4610.7660.776608.050856.133331.864Mianwali0.4160.7340.8230.7350.5030.3080.9180.7480.6250.3230.5880.679675.666978.739576.477Multan0.4890.7820.7950.6130.2920.3190.8720.6160.4250.3820.6370.674592.744365.6381005.901Muzaffarabad0.5370.7960.8370.5710.4250.1100.8910.7220.6000.4120.6360.655660.512348.687834.538Sialkot0.5230.7740.8220.7550.5920.4040.9230.7990.7080.4120.6270.662436.879238.422380.266SPTI-3Astore0.3050.5060.5580.9220.8720.8470.9780.9390.9250.2180.3770.425270.129115.508396.597Chhor0.2750.4310.4740.9050.8340.7810.9700.9230.9050.1900.3050.356266.512126.178351.577Chitral0.2520.4030.4610.9500.9170.9020.9850.9610.9510.1770.2960.34589.23646.108128.119Kalat0.3420.5040.5430.9020.8500.8240.9670.9270.9140.2350.3620.403202.44693.793164.865Kohat0.3950.5800.6250.8740.8110.7090.9590.7100.8800.2910.4310.473206.435125.313242.491Mianwali0.3300.5390.6320.9080.8330.7720.9720.9210.8890.2380.3940.472180.31976.808163.529Multan0.4160.6600.6840.8630.74606390.9520.8770.8430.2950.4850.522196.75092.480175.258Muzaffarabad0.3480.5220.5570.8250.8200.5220.9640.9170.8940.2630.3880.428249.060189.263490.751Sialkot0.3170.4760.5310.76307560.8010.9700.9270.9070.2290.3550.408230.092557.723263.315SPTI-6Astore0.2870.3510.4550.9450.9120.8940.9910.9810.9490.1980.2480.34775.91886.424129.332Chhor0.2550.3530.4300.9390.9160.8590.9790.9600.9360.1550.2290.30798.630131.632181.744Chitral0.2470.2990.4020.9640.9550.9250.9860.9780.9620.1670.2110.30681.65152.158131.476Kalat0.3160.3800.4390.9560.9420.8970.9930.9760.9520.1930.2400.30179.513300.351238.315Kohat0.8830.4550.5070.9220.9570.8230.9680.9830.9300.6950.3130.370173.514253.266703.016Mianwali0.3190.4090.4880.9120.8980.8660.9720.9530.9330.2180.2720.332209.50498.576422.293Multan0.3820.4780.4810.9220.8720.8410.9930.9630.9370.2640.3180.337179.67699.793136.787Muzaffarabad0.3900.4710.5270.8850.7530.3430.9590.9440.9140.2780.3370.395158.979114.439364.869Sialkot0.2980.3750.4320.8890.8820.6910.9750.9590.9440.1990.2490.314117.07897.213149.507SPTI-9Astore0.2880.3460.4340.9490.9360.8740.9800.9710.9510.1970.2480.30976.93362.029194.755Chhor0.2830.3550.4330.9430.9100.7640.9770.9510.9070.1580.2210.289150.82892.464129.091Chitral0.2630.3080.4430.9580.9550.9140.9840.9790.9570.1850.2220.319116.52758.64586.048Kalat0.2840.3440.4090.9510.9450.9270.9820.9740.9630.1690.2100.241115.838171.454201.090Kohat0.3180.3710.4150.9420.7650.8460.9750.8750.9390.2200.2590.30289.845123.001117.730Mianwali0.2640.3070.3880.9480.9480.9160.9820.9750.9610.1840.2190.277466.11896.627164.084Multan0.3660.4360.5470.9280.9050.8620.9680.9530.9280.2470.2940.352143.63490.969124.937Muzaffarabad0.3230.3850.4710.9380.9210.8490.9750.9640.9410.2400.2890.357180.985107.083712.336Sialkot0.2980.3340.4230.9470.9380.8950.9760.9700.9510.2020.2340.310121.5391223.05552.815SPTI-12Astore0.3110.3360.3000.9510.9450.9380.9730.9770.9780.2030.2160.204115.137158.271131.494Chhor0.3020.3400.3020.9460.9330.8640.9740.9670.9730.1510.1540.15897.432241.04060.678Chitral0.3480.3720.3400.9400.9360.9230.9710.9680.9730.2140.2180.21394.94283.428144.242Kalat0.2400.2840.2570.9730.9680.9660.9890.9850.9870.1220.1310.12451.69059.213457.499Kohat0.3010.3420.3020.9520.9440.9390.9780.9720.9780.2030.2240.204146.08178.430131.950Mianwali0.2540.2770.2540.9680.9590.3280.9850.9820.9840.1690.1800.16797.51894.95070.970Multan0.3420.3830.3620.9350.8140.8990.9570.8130.9620.2160.2290.225412.64480.412142.043Muzaffarabad0.3440.3700.3350.9410.9350.9270.9720.9680.9720.2460.2600.242111.26090.649165.662Sialkot0.3290.3640.3290.9380.9320.9110.9720.9660.9710.2090.2240.213108.491114.55492.435Similarly, the numeric quantities of WI at the Astore station are 0.999, 0.748, and 0.664, respectively. The KGE index indorses ELM’s superior performance at all the stations by providing maximum values as compared to MLP and ARIMA. Similarly, Willmott’s agreement “WI” index consistently provides the highest values for the ELM algorithm. KGE and WI are considered the most appropriate metrics for the performance assessment of meteorological and hydrological models. The KGE is calculated using the correlation coefficient, the ratio of variations, and the ratio of averages of predicted and observed series using equation (9). The values of “WI” for the ELM model for all the selected stations are close to 1, which endorses the ELM as the best performing model. The similar superior performance of ELM continued for other time scales at selected stations. Overall results for the training phase show that the ELM model has shown good agreement at all selected stations.As the time scale of the drought index increases, the performance of the proposed algorithm improves. As a result, the ELM algorithm showed superior performance to its competitive algorithms (MLP and ARIMA). A comparison of all the performance assessment metrics concluded that ELM algorithm is selected as the adequate model for the estimation and forecasting of drought indices (see Table7).The consistency and co-movement of the observed and estimated values of SPTI are further assessed while employing Karl Pearson’s product-moment correlation coefficient. Table8 shows the numerical results of the correlation coefficient between the observed and predicted values of SPTI using ELM, MLP, and ARIMA models for training data. The quantitative results of the Astore station for ELM, MLP, and ARIMA are 0.87, 0.59, and 0.56, respectively, indicating a better agreement of ELM to predict SPTI-1 contrary to other candidate algorithms. At Chitral station, values of correlation are (0.93, 0.74, and 0.69). For SPTI-3 at Astore station, these results are 0.96, 0.88, and 0.86, respectively, and at Chitral station, the correlation values for ELM, MLP, and ARIMA are 0.97, 0.93, and 0.90. These numerical results clearly illustrated that the performance of the ELM model significantly improved as the time scale increased. A similar pattern of the superior prediction performance of ELM continued for other time scales at all the selected stations. Although all the models have shown reasonable prediction performance, quantitative results evidently confirmed that the estimated values of SPTI using the ELM model strongly correlate with the observed values of SPTI for any time scale.Table 8
Pearson’s correlation coefficient between the observed and predicted values of SPTI 1–12 month time scales using ELM, MLP, and ARIMA algorithms at selected meteorological stations.
ScaleAstoreChhorChitralKalatKohatMianwaliMultanMuzaffarabadSialkotELMMLPARIMAELMMLPARIMAELMMLPARIMAELMMLPARIMAELMMLPARIMAELMMLPARIMAELMMLPARIMAELMMLPARIMAELMMLPARIMASPTI-10.870.590.560.860.620.540.930.740.690.860.630.520.840.480.390.880.590.430.840.410.320.840.560.480.880.660.60SPTI-30.960.880.860.950.860.830.970.930.900.940.860.840.920.830.800.950.850.800.910.780.750.930.850.820.940.860.83SPTI-60.960.940.900.960.920.880.970.960.930.950.930.910.940.900.870.940.910.870.930.890.890.930.890.860.950.920.89SPTI-90.960.940.910.960.280.840.970.960.920.970.950.930.950.770.890.960.950.920.940.910.860.950.930.890.950.940.91SPTI-120.960.950.960.950.930.950.940.940.950.980.970.980.960.950.960.970.960.970.770.920.930.950.940.950.950.930.94Usually, climatic and meteorological studies comprise high-frequency datasets that require fast algorithms. So, speed is a notable characteristic for determining the reliability of the algorithm. The algorithm selection for climatic studies is subjective in terms of speed and relative efficiency. ELM has the novelty of being the fastest algorithm among the ANN class to solve complex datasets. ELM algorithm training and testing time were almost 32 times faster than ANN, indicating ELM's supremacy over other ANN algorithms [58].The functional relationship between the actual (observed) and predicted values of SPTI using ELM and other algorithms for the “Astore” station is shown in Figure4 using a line graph, which depicts significantly less variation among the observed values of SPTI and the predicted values using ELM.Figure 4
Performance assessment of ELM, MLP, and ARIMA models through co-movement of observed and predicted values of SPTI at Astore station for the training phase.MLP and ARIMA were unable to capture all the shocks in historical values of SPTI, and departure from observed values was significant. Here, the ELM algorithm reflects more precise and accurate predictions. Figures5– detail the ELM's prediction performance for all the time scales at the selected stations, which depicts the ELM algorithm's better prediction performance.Figure 5
Co-movement of the observed and predicted values of SPTI using ELM at Astore, Chhor, and Chitral stations for the training phase.Figure 6
Co-movement of the observed and predicted values of SPTI using ELM at Kalat, Kohat, and Mianwali stations for the training phase.Figure 7
Co-movement of the observed and predicted values of SPTI using ELM at Muzaffarabad, Multan, and Sialkot stations for the training phase.It was observed that as the time scale increases, the prediction performance of algorithms substantially improves. For SPTI-12, significantly fewer deviations have been observed among drought index’s paired (observed and predicted) values. These multi-line plots indicate that the ELM model incurred smaller errors than the two counterparts. The predicted values using ELM model follow the observed values of SPTI more precisely. These graphical presentations show that the ELM model has shown more accuracy than MLP and ARIMA models for estimating SPTI for various lead time scales at selected meteorological stations.Scatter plots of the observed and predicted time series data sets are another way to assess the prediction performance of probabilistic, machine learning, and ANN algorithms. Figures8–10 show the scatter plots of the historical observed and predicted values of SPTI-6 using ELM, MLP, and ARIMA models at all the selected meteorological stations. Predicted values through the ELM model have shown a strong correlation with the observed values of SPTI through scatter plots. Another graphical presentation endorses the superior performance of ELM. We can observe that the ELM algorithm showed more accuracy and can potentially predict drought conditions in any climatic zone.Figure 8
Scatter plots of the observed and predicted values of SPTI-6 using ELM, MLP, and ARIMA at Astore, Chhor, and Chitral stations for the training phase.Figure 9
Scatter plots of the observed and predicted values of SPTI using ELM, MLP, and ARIMA at Kalat, Kohat, and Mianwali stations for the training phase.Figure 10
Scatter plots of the observed and predicted values of SPTI using ELM, MLP, and ARIMA at Muzaffarabad, Multan, and Sialkot stations for the training phase.Figure11 represents the Taylor diagrams for the Astore station with (1, 3, 6, 9, and 12) month time scales for the training phase. Taylor diagrams are the more comprehensive and precise way to represent the estimation and forecast ability of a model. The similarity between the predicted and actual values of SPTI is evaluated in terms of correlation (as a measure of closeness), and the variation is assessed by the standard deviation (SD) and the RMSE. For SPTI-1, the correlation of modeled data using the ELM algorithm with actual observations was about 0.9, followed by MLP and ARIMA with 0.6 each. As time scale increases, the prediction performance of algorithms significantly improves. For SPTI-3, the ELM algorithm is significantly closer to the actual values as its correlation is about (0.97) as compared to MLP (0.9) and ARIMA (0.85). The Taylor diagram exhibits the superior performance of ELM algorithm for estimating SPTI with (1, 3, 6, 9, and 12) month lead time scales (see Figure 11). Figure 12 illustrates the violin plots related to the training phase of ELM and other candidate models. The red dot represents the mean value, a thick white bar represents the interquartile, and a thin blue line represents the whole data set distribution. These are the components of the boxplot, but the colored area on both sides of the blue line is the estimated kernel density of the shape of the distribution of data. Higher probabilities are associated with wider parts of the diagram, and thinner sections show lower probabilities. This is another systematic way to compare the prediction performance of models. It is noticeable that the mean of the predicted and observed values of SPTI-6 for the ELM model was almost similar for all the stations, whereas other models have slight variations. All the models have shown reasonable agreement but Violin diagrams of ELM and observed data are nearly identical. These graphical illustrations affirmed that the proposed ELM model is better in estimating the actual values of SPTI at all selected meteorological stations.Figure 11
Taylor diagrams for the evaluation of the prediction performance of ELM, MLP, and ARIMA algorithms at Astore station with different time scales in the training phase.Figure 12
Violin plots of the observed SPTI-6 and the modeled data using ELM, MLP, and ARIMA models at all the selected meteorological stations in the training phase.After calibrating and validating algorithms for training datasets, the proposed algorithm's generalization capability has been assessed in the next step. Finally, out of the sample forecast of SPTI for all the lead time scales, (1–12) is carried out for 36 months from 2014 to 2016. These forecasts are considered sufficient for drought preparedness and mitigation policies. Similar performance metrics have been used to analyze the difference between the observed and forecasted values of SPTI for different time scales at all selected meteorological stations. Numerical results related to these performance metrics are given in Table9. If the comparison of KGE was made at Astore station, the values of KGE for ELM, MLP, and ARIMA models are 0.988, 0.771, and 0.671, respectively. The KGE for the ELM model is significantly higher than its counterparts for all the time scales at all the selected meteorological stations, which clearly endorses the better forecast capability of the ELM model. Similarly, WI and MAE quantitative results also ratified that the ELM model outperformed the ANN and ARIMA models. RMSE and MAPE endorsed ELM as a better forecasting model for most stations.Table 9
Forecast performance of ELM, MLP, and ARIMA algorithms for SPTI with 1–12 month time scales at nine meteorological stations.
StationRMSEKGEWIMAEMAPE (%age)ELMMLPARIMAELMMLPARIMAELMMLPARIMAELMMLPARIMAELMMLPARIMASPTI-1Astore0.6881.0340.8930.9880.7710.6710.9970.8970.8840.7980.8550.865247.353251.284364.431Chhor0.5790.6000.5590.8710.7120.5410.9990.7300.6630.5310.4470.468559.287269.288396.060Chitral0.6430.6470.6900.9090.6850.6780.9790.8360.8230.4950.5100.520142.102193.952414.532Kalat0.7760.8320.7270.9520.7790.7670.9980.9660.9760.6760.7310.684500.369824.503616.762Kohat0.9560.9701.0310.9790.9300.9000.9990.9920.9540.7910.8150.869528.0661136.164268.277Mianwali0.7730.9840.8010.9610.7530.6900.9980.9540.9230.6520.8180.655309.146369.444860.231Multan0.7140.6660.8240.9860.4280.6870.9970.7400.8850.6120.5770.706545.5681291.2043958.975Muzaffarabad0.7640.7430.7240.9370.8770.6730.9970.9660.9940.6050.5970.553321.8103469.574455.340Sialkot0.7200.7450.8220.9310.9020.9140.9950.9890.9820.5540.5930.66391.037151.493269.868SPTI-3Astore0.7640.8960.6070.9700.7440.5110.9910.8650.7600.5890.7350.488797.80063.99284.646Chhor0.5630.5030.5390.86609620.6780.8210.8730.8710.4380.3660.385193.072159.359141.894Chitral0.4680.5220.5140.9210.7090.8200.9650.9080.9180.4080.3960.4052354.0679.394150.935Kalat0.9871.2471.1520.9960.8580.6780.9990.9560.7580.8501.0670.899329.382314.217199.070Kohat0.7710.7750.8730.9740.9650.730.9940.9930.8910.7290.6680.726116.387330.497134.998Mianwali1.1041.2211.0410.9970.9610.9430.9970.9870.9560.9141.0140.8061082.42561.266173.703Multan0.8010.7900.9820.8940.9270.9380.9990.9970.9780.6600.6540.826210.083193.947210.070Muzaffarabad0.7910.7850.7770.9670.8560.8730.9990.9900.9960.6020.6860.601230.935809.403334.650Sialkot0.6810.7250.7680.9380.9240.6280.9960.9770.8110.5200.5320.601265.276259.027201.413SPTI-6Astore0.4020.6700.4820.9620.6070.5410.9850.7820.7330.4920.5820.40846.301343.20737.464Chhor0.8190.8890.6690.9100.6870.4560.9980.9670.8760.7020.7710.475697.320269.54173.976Chitral0.3740.4050.3990.9610.8620.8460.9840.9350.9430.3080.3250.308172.828140.348140.661Kalat0.9720.7921.2310.9650.9210.9330.9930.8840.9760.8070.6661.0221364.65335.608184.772Kohat0.9570.9300.8690.9240.8760.7590.9900.9480.9180.7640.8410.71989.183301.444109.716Mianwali1.4112.1481.4570.9740.7880.8850.9990.9750.9871.3191.8681.145347.413215.620298.005Multan0.7831.2091.0760.9310.7810.9260.9990.9990.9910.6730.9820.886222.776194.657173.698Muzaffarabad0.7711.1760.8720.9890.9350.8930.9990.9950.9810.6480.8220.6952458.04214.1594533.34Sialkot0.8800.7330.6560.9530.9440.9280.9980.9910.9940.7560.6700.557776.441532.295178.914SPTI-9Astore0.6340.6470.8740.9680.7280.7830.9900.8570.8690.4960.5420.77641.23049.977101.519Chhor1.1341.3280.9610.9490.9360.9450.9940.9960.8290.9441.1510.760358.134243.020584.953Chitral0.3720.5040.3830.9530.8630.8200.9840.9210.9100.3220.4060.29753.6963157.8534.444Kalat0.4281.1330.7820.9910.7240.8450.9990.8920.9350.8760.9820.713187.052360.12550.731Kohat1.0321.0420.8360.9410.9320.7740.9920.9870.9620.8400.7910.61798.19087.81372.893Mianwali1.3502.0621.4570.9700.9240.6160.9980.9820.9641.4291.8951.145253.908162.066298.005Multan0.7460.6901.1480.9730.8730.8980.9990.9980.9990.6090.5950.9461139.221061.42134.464Muzaffarabad0.7630.8000.7610.9850.9480.9380.9990.9790.9900.6340.6600.59278.798144.12688.432Sialkot0.7650.6810.8390.9250.9040.8510.9970.9960.9850.5890.4850.65594.806227.032221.098SPTI-12Astore0.2580.7541.2290.9620.7210.7700.9920.9060.8690.2150.6581.20210.93668.122187.249Chhor0.6700.7070.7280.9870.9750.8790.9960.9350.9210.5660.6000.648197.9181307.502828.820Chitral0.5210.5240.7720.9960.7470.6780.9990.9220.8540.4060.4130.65338.71035.631117.257Kalat1.3461.2161.1680.9750.9940.6120.9980.9980.8041.1791.0260.9841739.84558.414343.401Kohat0.6151.0660.9630.9810.9590.9000.9970.9950.9810.8600.7900.81984.52266.743144.398Mianwali1.4021.5131.4840.9940.8820.7750.9990.9890.8371.6291.2981.273147.015180.588169.143Multan0.5611.0720.6640.9950.9160.8700.9980.9740.9700.4340.9970.55360.9761233.10223.381Muzaffarabad0.7201.1170.7540.9850.7560.8900.9970.8990.9620.5640.7800.58649.482257.15694.969Sialkot0.7930.7270.7660.9450.9260.9220.9930.9920.9730.6300.5930.59361.040181.0031734.785The functional relationship between the observed and forecasted values of SPTI using ELM, MLP, and ARIMA models for all the time scales at Astore stations is illustrated in Figure13 for the test phase starting from January 2014 to December 2016 for 36 months. A significantly smaller degree of deviation in SPTI for the ELM model and observed values of SPTI were exhibited. Although MLP has shown a reasonably good forecast performance compared to the stochastic seasonal ARIMA model, the ELM model evidently shows superior forecast performance.Figure 13
Observed and forecasted SPTI by ELM, MLP, and ARIMA models at Astore station for 1, 3, 6, 9, and 12 month lead time scales.In order to check the appropriateness of the ELM model, scatter plots were prepared using time series data of the observed and forecasted SPTI at Astore station for all lead time scales (Figure14). These scatter plots have shown a significant difference in the forecast performance of ELM and other models. The scatter plots depict the correlation, goodness-of-fit, and the extent of agreement between the observed and forecasted SPTI. The ELM model also clearly outperformed MLP and ARIMA models for the testing phase for all the selected time scales.Figure 14
Scatter plots of the observed and forecasted values of SPTI using ELM, MLP, and ARIMA models at Astore station with all lead time scales for the test phase. The least-squares regression line is also included.


## 3.3. Discussion

Drought is a multifaceted and commonly occurring hazard in several parts of the world. Its impacts are prevalent in the agriculture, socio-economic, and energy sectors. However, precise drought monitoring and estimation techniques can assist in decreasing the vulnerability of society to drought. The primary objective of the current study was to test the appropriateness and usefulness of the ELM model relative to other ANN (MLP) models and stochastic (ARIMA) models for drought forecasting. The prediction and forecast performance of ELM is compared with other ANN algorithms (MLP) and statistical stochastic (ARIMA) models. The prediction and forecast performance of models is assessed using numerous performance metrics, including RMSE, MAE, MAPE, KGE, WI, and Karl Pearson's correlation coefficient. The quantitative assessment revealed that both the ANN models (ELM and MLP) performed better than the stochastic model (ARIMA), and among the ANN models, ELM has shown supremacy by producing the smallest RMSE, MAE, and MAPE values and the maximum values for KGE, WI, and correlation coefficient for almost all the meteorological stations.Furthermore, ELM shows better agreement for both the training and test phases to predict the SPTI at all climatic stations than its counterparts. The efficiency of ELM, contrary to other models, is evident based on the performance metrics. A similar forecast performance has continued for higher-order time scales, consistent with earlier studies [30, 58]. Computational time consumed by drought modeling algorithms also needs to be optimized. Usually, large datasets are used as input variables for real-time drought modeling, which affects the computational performance of different models in terms of time. The ELM model is significantly faster than its counterparts. By evaluating all the numerical results of performance metrics and different graphical illustrations, it can be easily concluded that the ELM model attains the most accurate drought forecasting performance during training and test phases. The study revealed that ELM is the most appropriate, reliable, and efficient algorithm for drought prediction and forecasting. This study suggests that ELM can be used as an early warning drought forecasting tool for developing drought mitigation policies.


## 4. Conclusions

The reliable, efficient, and faster drought forecasting algorithms are useful for freshwater resource managers and drought mitigation policymakers. The current study examines the forecast performance of new machine learning (ELM) model using the Standardized Precipitating Temperature Index (SPTI). For application, meteorological time series datasets of monthly precipitation and minimum and maximum temperatures were collected from nine meteorological stations located in various climatological zones of Pakistan. Further, the prediction and forecast performance of the ELM model was compared with MLP and ARIMA models using different statistical performance metrics and graphical illustrations. The primary objective of the study was to investigate the appropriateness of the ELM model for predicting the nonlinear and complex temporal behavior of SPTI. The ELM model outperformed MLP and ARIMA models by producing the smallest root mean square error, mean absolute error, mean absolute percent error values, and maximum values for KGE, WI, and Karl Pearson’s correlation coefficient for all the selected meteorological stations for different selected time scales. KGE and WI unanimously endorsed ELM as the best forecasting model at all the stations. Moreover, by comparing forecasting results for a one-month time scale, RMSE clearly affirmed ELM as superior forecasting model for five stations including Astore (RMSE = 0.688), Chitral (RMSE = 0.643), Kohat (RMSE = 0.956), Mianwali (RMSE = 0.773), and Sialkot (RMSE = 0.720). Furthermore, MLP better performed at Multan (RMSE = 0.666) and ARIMA at three of the meteorological stations containing Chhor (RMSE = 0.559), Kalat (RMSE = 0.727), and Muzaffarabad (RMSE = 0.724). While using MAE as a performance measure, the ELM algorithm performs better for six [94] stations (Astore 0.798, Chitral 0.495, Kalat 0.676, Kohat 0.791, Mianwali 0.652, and Sialkot 0.554), MLP for two [95] stations (Chhor 0.447 and Multan 0.577), and ARIMA for one of the stations (Astore 0.765). MAPE recognizes [96] ELM as an appropriate algorithm [97] for seven stations (Astore 247.3, Chitral 142.1, Kalat 500.3, Mianwali 309.1, Multan 545.5, Muzaffarabad 321.8, and Sialkot 91.0) and MLP (Chhor 269.2) and ARIMA (Kohat 268.2) for one station. Contrary [98] to MLP and ARIMA models, the ELM model has [99] the super-fast computation capability of [100] drought modeling. These performance [101] comparisons clearly ratified the novelty and [102] appropriateness of the proposed ELM algorithm.In summary, this study suggests that ELM can be used as an early warning drought forecasting tool for developing drought mitigation policies using time series data of the Standardized Precipitating Temperature Index. The scope of the study can be enhanced by using wavelet data mining transformation to get the optimized forecast performance of the ELM model. ML algorithms have certain limitations, like requiring fast computing hardware. ANN models (ELM and MLP) have complex network structures among the machine learning algorithms. Another limitation of ANN algorithms is that there is no specific rule to determine the final network structure. Instead, the appropriate network structure is finalized through trial-and-error-based.



---

*Source: 4998200-2022-09-20.xml*