The Development of a Hybrid Wavelet-ARIMA-LSTM Model for Precipitation Amounts and Drought Analysis

: Investigation of quantitative predictions of precipitation amounts and forecasts of drought events are conducive to facilitating early drought warnings. However, there has been limited research into or modern statistical analyses of precipitation and drought over Northeast China, one of the most important grain production regions. Therefore, a case study at three meteorological sites which represent three different climate types was explored, and we used time series analysis of monthly precipitation and the grey theory methods for annual precipitation during 1967–2017. Wavelet transformation (WT), autoregressive integrated moving average (ARIMA) and long short-term memory (LSTM) methods were utilized to depict the time series, and a new hybrid model wavelet-ARIMA-LSTM (W-AL) of monthly precipitation time series was developed. In addition, GM (1, 1) and DGM (1, 1) of the China Z-Index (CZI) based on annual precipitation were introduced to forecast drought events, because grey system theory specializes in a small sample and results in poor information. The results revealed that (1) W-AL exhibited higher prediction accuracy in monthly precipitation forecasting than ARIMA and LSTM; (2) CZI values calculated through annual precipitation suggested that more slight drought events occurred in Changchun while moderate drought occurred more frequently in Linjiang and Qian Gorlos; (3) GM (1, 1) performed better than DGM (1, 1) in drought event forecasting.


Introduction
As one of the most destructive natural calamities, drought occurs when rainfall amounts are below normal for a long period. The characteristics are high frequency, long duration, wide influence [1,2], and damaging effects on grain yields and water supplies, so it is of great significance to model and forecast the rainfall amount and drought. Accurate precipitation predictions are required for the precise estimation of drought in an area [3]. More accurate and timely rainfall prediction can boost drought research, while greatly improving future water management policies in many ways.
Due to the nonlinear, stochastic and highly complex nature of rainfall data, timely and exact rainfall forecasting has remained a challenging task, and more complex technologies are needed. The autoregressive integrated moving average (ARIMA) and neural network (NN) are broadly trending [1,4,5]. ARIMA has good prediction accuracy and flexibility for different types of time series data, such as those found in hydrology [6,7], finance [8], agriculture [9], and medicine [10]; however, ARIMA cannot adequately simulate 1) has the advantages of fully fitting pure exponential sequences, no restrictions on the development coefficient, and broadens the application scope of the model. However, in actual situations, the unique superiorities of DGM (1, 1) cannot be fully utilized because the data are basically inconsistent with exponential growth, which causes scholars to usually choose GM (1, 1) instead of DGM (1, 1) when solving practical problems. Since first being proposed, grey prediction models have been broadly employed in various fields, such as electricity consumption prediction [40,41], air pollution forecasting [42][43][44] and energy forecasting [45][46][47]. Such practical applications show that grey prediction models have wide applicability, especially in situations with incomplete information and inaccurate data. Grey prediction models have successfully dealt with various problems, but only a few scholars have used these models to study drought prediction, while prediction of drought events conforms to the characteristics of grey systems. Therefore, this paper uses GM (1, 1) and DGM (1,1) to predict the occurrence of drought events. In addition, we choose GM (1, 1) and DGM (1,1) because that they are the most basic and widely used grey prediction models.
Generally, hybrid models have higher prediction accuracy than single models. In this paper, the wavelet-ARIMA-LSTM method are proposed for the first time. It combines the advantages of wavelet, ARIMA and LSTM and can predict future precipitation more accurately. Drought is considered to be the most incomprehensible and least understood disaster by many researchers. It is very difficult to achieve accurate forecasting of drought events when facing the problems of insufficient samples and poor information. The grey system model places a particular emphasis on dealing with the uncertainty brought by small samples. Based on this, the study uses GM (1, 1) and DGM (1, 1) to predict drought years for drought risk analyses and drought warnings. The major objectives of this study are: (1) to develop a hybrid wavelet-ARIMA-LSTM method to predict monthly precipitation for the period 1967-2017 in Northeast China; (2) to analyze drought characteristics in Northeast China based on the drought index, CZI; and (3) to use GM (1, 1) and DGM (1,1) to predict the occurrence of drought events and compare the predictive capabilities of the two methods. It is expected that the research results will help to provide decision support for rainfall predictions, which in turn will help in planning adaptative measures to reduce drought impacts and provide decision support for disaster prevention.

Study Area
In this study, the study area includes three stations, namely Changchun, Linjiang and Qian Gorlos in Jilin province, Northeast China between 41 • N to 46 • N and 122 • E to 131 • N and experience a temperate continental monsoon climate. The climate of Jilin province is classified into different categories of humid climatic conditions. Qian Gorlos is located in northwestern Jilin province, which is arid and semi-arid, while Linjiang is located in southeastern Jilin province, which is humid and semi-humid. The climate of Changchun represents a transitional zone between the semi-humid mountains to the east and semi-arid plains to the west. The locations of the stations used in Jilin province, Northeast China, shown in Figure 1. Table 1 presents the geographical coordinates and climatic conditions of the three selected stations.

Site Precipitation Observations
The observed monthly and annual precipitation data for the studied regions were collected from the National Meteorological Information Center (NMIC) of the China Meteorological Administration (CMA) from January 1967 to December 2017. In this study, monthly rainfall time series data from the studied stations were utilized for precipitation predictions, while annual rainfall amounts were utilized for drought analyses and predictions. For precipitation predictions, the monthly data between 1967 and 1997 (approximately 60% of the total data, i.e., 31 × 12 = 372 data points) were used to train the models, and the monthly data from 1998 to 2017 (40% of the total data, i.e., 20 × 12 = 240 data points) were used to test the models. For drought analyses and predictions, the annual data between 1967 and 1997 (approximately 60% of the total data, i.e., 31 data points) were employed to train the grey prediction models, and the remaining data were employed to test these models. Figure 2 shows the time series plots of observed monthly precipitation over 51 years for three stations throughout the study period. As indicated in Figure 3, precipitation in July and August is much higher than in other months, and the average value for each month is not the same. These data reflect a notable seasonal effect, which is consistent with the sequence chart shown in Figure 2. It is worth noting that precipitation time series must be normalized to eliminate the dimensions of observational precipitation datasets and map the data to the range of 0~1, which is more convenient and faster. Here, all monthly precipitation data were normalized as follows: where x , x, x min and x max denote the normalized precipitation data, observed precipitation, minimum value of observed data, and maximum value of observed data, respectively.   The ARIMA method proposed by Box and Jenkins [48] has gained great popularity in many fields, and research experience has confirmed its strength and flexibility [1,10]. It is a stochastic sequential model that is trained to forecast future data points. The model can capture complex patterns and relationships as it can combine capturing observations of lagged terms and white noise. The ARIMA consists of three parts: autoregressive (AR), integration (I), moving average (MA). The corresponding parameters are p, d, and q. The general ARIMA model is called ARIMA (p, d, q). The method is composed of three main steps: identification, estimation parameters, and forecasting [1,4,49].

Long Short-Term Memory Method (LSTM)
A traditional RNN which is a type of artificial neural network, readily introduces the problems of gradient disappearance and explosion, thereby making it difficult to capture long-term time correlations. Long short-term memory (LSTM) is a type of time-cyclic neural network, that is specifically used to solve the long-term correlation problem of general RNN. The LSTM proposed by Hochreiter and Schmidhuber [50] was initially used in the field of deep learning and was popularized by researchers in subsequent work [51].
The network takes three inputs and two outputs, as shown in Figure 4. For the inputs, x t is the input of the current time step, h t−1 is the output of the last LSTM unit, c t−1 is the memory of the previous unit, h t is the output of the current network, and c t is the memory of the current unit. The LSTM model has an input gate i t , output gate o t and forget gate f t . There are three stages in the LSTM. The first is the forgetting stage, which mainly is used to selectively forget the input from the previous node. Specifically, the calculated f t is used as the forget gate to control which parts of c t−1 in the previous state need to be retained and which needs to be forgotten. The second stage is the selective memory stage, which selectively memorizes the input x t . If input x t is important, it should be noted down, and if it is not, it should be noted less. The third stage is the output phase, which determines which outputs will be treated as the current state [52]. The LSTM equations are as follows [22,52]: where f t , i t and o t present the activations of the forget gate state, input state and output gate, respectively, at time step t; c t is the current input cell state; c t and c t−1 are the cell state vectors at time t and t − 1; h t and h t−1 are the hidden state vectors also known as output vectors at time t and t−1; σ and tanh denote the sigmoid function and hyperbolic tangent function; W f and b f represent the weight matrix and bias of the forget gate layer; similarly, W i and b i represent the weight matrix and bias of the input gate, W c and b c represent the weight matrix and bias of the unit state, and W o and b o represent the weight matrix and bias of the output gate, respectively.

Discrete Wavelet Transformation (DWT)
Wavelet analysis methods can be classified into continuous wavelet transformation (CWT) and discrete wavelet transformation (DWT) [1,53]. The main difference between the two is that continuous transformations operate on all possible scaling and translation values, while discrete transformation uses a specific subset of all scaling and shifting values. The main disadvantage of the CWT is that the construction of the CWT inverse is more complicated and thereby computationally difficult. Discrete wavelet transforms are widely used in the prediction field because of their short calculation times and easy application. This paper chooses DWT since DWT simplifies the transformation process and reduces the workload; the discrete wavelet transform can still produce very effective and accurate analysis results. DWT adopts the following form [53]: where a and b are integers that control the scale and time; ψ denotes the mother wavelet; s 0 denotes a dilation step with a constant value that is greater than 1; and γ 0 represents a position variable with value greater than zero. The most common selections for the parameters for s 0 and γ 0 are 2 and 1, respectively. When a time series is discrete with a value of x t occurring at discrete time t, the wavelet coefficient (W Ψ (a, b)) of DWT becomes [53]: The wavelet coefficients of the wavelet transform are calculated at scale s = 2 a and locations γ = 2 a b, which reveal the signal changes at different scales and locations [53].

Development of Wavelet-ARIMA-LSTM (W-AL) Model
In time-series applications, although there are many available time series models, none of them can provide the best results in various situations. A large number of time series prediction studies have indicated that hybrid methods can improve prediction performance [4]. By making full use of the advantages of each method in the combination model, the error risk from using an inappropriate method is reduced, and more accurate results are obtained. In this study, we develop a new hybrid method for time series forecasting that combines the strengths of wavelet transformation, ARIMA and LSTM.
The method is divided into decomposition and reconstruction. First, the original sequence is decomposed by high pass (detail) and low (approximate) pass filters, and the high-frequency and low-frequency components of the sequence are extracted respectively. Then, ARIMA is used to estimate the approximate signal, and LSTM is used to estimate the detailed part of the signal. Finally, the predicted wavelet coefficients obtained are used to reconstruct the data. The main advantage of WT is that it can analyze and process over different time scales. Figure 5 shows the framework of the wavelet-ARIMA-LSTM model development. Debauches' (db4) mother wavelet was used for decomposing the rainfall time series in this study. The monthly rainfall prediction accuracy of the W-AL models was compared to that of the single ARIMA and LSTM models.

Evaluation Metrics
To analyze the reliability and forecasting performance of the model, it is necessary to verify the accuracy of the models. The root mean square error (RMSE), which represents the standard deviation of the predicted results of the models; mean absolute error (MAE) which directly provides the average difference between the predicted results and actual data; and coefficient of determination (R 2 ), which provides a way to assess the results of the same model on different data, are adopted to assess the performance of the models. Smaller RMSE and MAE values indicate better model performance, and larger R 2 values indicate better model performance. The criteria are defined as follows: where n, x(i),x(i) and x represent the number of observations, observed data, predicted data and the mean of the observed data, respectively.

Grey System Models on Drought Events
There are many kinds of indicators for assessing drought, and the annual precipitation amount is an important sign of drought. Drought determined by using annual precipitation is generally called meteorological drought, and the year in which the meteorological drought occurs is referred to as the drought year. This paper chose the drought index CZI.
CZI, which was developed by the National Climate Centre (NCC) of China in 1995 as an alternative to the SPI [30] is used to describe drought conditions [28,29,31]. The value of CZI is calculated as: where C s is the coefficient of skewness, ϕ i is the standard variation, and the calculation formulas can be represented as follows: where x is the mean of the observation values, and n is the number of observation values. The classifications of drought severity levels for CZI are given in Table 2.

Grey Prediction Models
Grey prediction theory, as a significant part of grey system theory, addresses problems with small sample sizes and inadequate information [37]. The most important feature of grey prediction models is that they have relatively loose requirements for the collected data for the study of the problem. This theory can take all random variables as the object of study, and then regard their random nature as a time-related grey process. When grey prediction models are applied to the prediction of drought years, there is no need to know any a priori characteristics of the original data distribution, the test accuracy after modeling is high, and the models can better reflect the actual situation.
Drought is a complicated phenomenon and is one of the most unpredictable natural disasters. Therefore, data collection and selecting the potential influencing factors of drought are difficult tasks. Meanwhile, drought occurrences are irregular and discontinuous events, and their prediction methods are more difficult. Grey prediction models can use fewer data to obtain the desired results. Thus, considering that predictions of drought years conform to the characteristics of grey system models, this paper used GM (1, 1) and DGM (1, 1) [38,39,47], as the most fundamental and extensively used grey prediction models. The descriptions of the processes and computations of the GM (1, 1) and DGM (1, 1) methods are detailed in Wang et al. [38].
GM (1, 1) has the problem of prediction error caused by the abrupt change from discrete to continuous. Therefore, DGM (1, 1) proposed by Xie and Liu [38] make up for the defects of the traditional GM (1, 1). DGM (1, 1), which can be called the discrete form of the GM (1, 1) model, is superior to the GM (1, 1) because it can fully fit the pure exponential sequences and has no limit to the development coefficients. However, due to the interference of random factors in the actual data generation process, the superiority of DGM (1, 1) over GM (1, 1) cannot be widely and reliably verified in practical applications. For the univariate non-negative time series x (0) = (x (0) (1), x (0) (2), . . . , x (0) (n)).
The sequence x (1) = (x (1) (1), x (1) (2), . . . , x (1) (n)) is the first-order accumulation generation of x (0) , where x (1) (k) = ∑ k i=1 x (0) (i), k = 1, 2, · · · , n. The DGM (1, 1) is defined as follows [38]: Tests of grey prediction models mainly include residual tests and posterior variance tests. The residual test calculates the absolute e(i) = x (0) (i) −x (0) (i), (i = 1, 2, · · · , n) where x (0) (i) andx (0) (i) represent the univariate nonnegative time series and the first-order accumulation generation of x (0) respectively, and the relative error ε(i) = e(i) × 100%, (i = 1, 2, · · · , n) between the original sequence and the grey prediction sequence. The smaller the relative error, the higher the accuracy of the model. The posterior-variance test includes two indices: the variance ratio C and small error possibility P. The specific functions can be expressed as [38]: Then, the variance ratio C = S 2 S 1 , is computed and then the small error probability P = p(|e(i) − e| < 0.6745S 1 ) is computed. The accuracy of the models is determined according to the Table 3. If both the residual test and posterior variance test are qualified, the model can be used for prediction.

Time Series Analysis of Monthly Precipitation Amounts
The rainfall time series data are separated into two detailed subseries and one approximate subseries by db4 mother wavelet, which is a frequently used wavelet for the DWT. Because the two detailed subsequences of the wavelet are nonlinear, LSTMs are used to predict the nonlinear time series. Meanwhile, because one approximate subsequence of wavelet is linear, ARIMA is used for prediction. Finally, the predicted values of ARIMA and LSTM are used to reconstruct the data series.
In the present research, the proposed W-AL was compared with the single ARIMA and LSTM models by adopting three statistical indicators for evaluating the performance of W-AL for predicting precipitation at the monthly scale and the results of the comparisons for the test stages are presented in Table 4. Performance comparisons of the best-fitted models indicated that W-AL was the best performer and was followed by LSTM and then ARIMA.  21.712, and R 2 from 0.366 to 0.670. These results indicated that the best RMSE and MAE values were found for Qian Gorlos, while the best R 2 value was found for Linjiang. However, RMSE and MAE values have limitations; that is, the same algorithm model that is utilized to predict monthly precipitation at different stations, cannot reflect the fitting effect of the model at different stations. Because the dimensions of the data are different at different stations, it is impossible to directly compare the predicted values, and it is impossible to determine the stations for which the model performs better. In contrast, R2 converts the predicted results into accuracies, and the results all fall between 0 and 1. For the prediction accuracies of different stations, it is possible to employ R 2 to compare and determine which stations perform better. Based on this, it can be found that the fitting effect of the model was best in Linjiang with humid and semi-humid climate type and was worst in Qian Gorlos with arid and semi-arid climate type which revealed the model predicted better in humid region and worse in arid region.    In order to illustrate the superiority of the W-AL model, this paper also used 70:30, 80:20 and 90:10 as different proportions of training and test sets. The RMSE values also show that the proposed W-AL model is superior to LSTM and ARIMA under different ratios as shown in Table 5. It can be found that the prediction accuracies of the LSTM and W-AL present an overall trend of increasing with the increase of the proportion of training sets. The prediction accuracy reaches the highest when the ratio training:test is 80:20 at Changchun and Linjiang, while it is 90:10 at Qian Gorlos. More accurate rainfall prediction can not only boost drought research, but also become an important reference for the impending severe weather warning. The study area, Northeast China, is one of the most important grain and animal husbandry production regions, while modern statistical analyses on precipitations and droughts here are relatively limited. Northeast China, which is recognized as a sensitive area of climatic change in global climate models due to its continental monsoon climate, suffers drought events and rainstorm consecutively [54,55]. This has resulted in numerous negative impacts on the national economy of the region. Therefore, the investigation on the rainfall prediction play a vital function in improving the risk management and prevention of meteorological disaster such as droughts. In this paper, the W-AL can improve the prediction accuracy of monthly precipitation compared with single ARIMA and LSTM methods, and becomes a new method for the statistical procedures to predict rainfall. On one hand, ARIMA has better prediction accuracy and flexibility for different types of time series data than other linear methods [6][7][8][9][10], while linear methods cannot capture the nonlinear characteristics of rainfall processes. LSTM, as a nonlinear method, can automatically learn complex time patterns through high-level abstraction and nonlinear transformation and achieve approximations of complex functions compared with the simple NN models. Thus, LSTM has a good accuracy, especially when the variables are nonlinear over other nonlinear methods like GEP which is more sensitive to the quality of the measured data [18][19][20]. On the other hand, hybrid methods can combine the superiorities of several models to overcome the shortcomings of each single model, and accordingly improve prediction accuracy in most cases [1,4,24]. After all, this paper uses the W-AL, ARIMA and LSTM to predict monthly precipitation in Northeast China to illustrate the predictive power of the W-AL over any single model at different climate types. Besides, this paper considers different ratios of training and test sets to conduct further research on the superiority of the W-AL model.

Identification of Drought Events
Drought events at the three stations were analyzed by using CZI calculated from the annual precipitation time series data. Figure 8 displays the annual CZI that was obtained for Changchun, Linjiang and Qian Gorlos for the period from 1967 to 2017. As seen in

Projections of Drought Events
The GM (1, 1) and DGM (1, 1) prediction models were established for the annual drought events of Changchun, Linjiang and Qian Gorlos from 1998 to 2017 to predict the annual drought events from 1998 to 2017 and to compare these with actual drought events. According to the classification results of drought statistics in CZI, drought years for a 31-year period were selected. Since the grey prediction models have better prediction effects for small numbers of sample data, the corresponding numbers of drought years after 1972 were selected to establish the initial sequence in Changchun, Linjiang and Qian Gorlos. Finally, the parameter estimations of GM (1, 1) and DGM (1, 1) prediction models are shown in Table 6. The actual drought years and forecasted drought years of the two models are presented in Table 7. As seen, the GM (1, 1) model predicted that the drought years after 1998 in , which were the same as the GM (1, 1) results. In contrast, GM (1, 1) performed better than DGM (1, 1). The average relative error values for GM (1, 1) were higher in Linjiang and were followed by Qian Gorlos, while Changchun showed a minimum. In summary, the prediction results of GM (1, 1) performed better than DGM (1, 1) such that GM (1, 1) is used to predict drought years for different drought levels. Furthermore, GM (1, 1) predicted poorly in relatively humid regions and well in relatively arid regions.
This paper studies the rainfall and drought in Northeast China from two aspects. Firstly, take the monthly rainfall time series data into consideration, and the proposed hybrid model W-AL is employed to predict the rainfall, so as to improve the prediction accuracy and carry out drought warning better, which has been discussed in the previous section. Secondly, considering the annual rainfall data, drought events and severity can be identified by drought index CZI and the occurrence of drought events can be predicted by grey system methods GM (1, 1) and DGM (1, 1). The grey prediction models have been widely concerned by academic circles [35][36][37][38][39][40][41][42], but only a few scholars have used the models to study drought prediction. The biggest characteristic of the grey prediction model is that the data collected are relatively loose, which solves the problems of small sample size and insufficient information. Drought occurrence, which is irregular and discontinuous, conforms to the characteristics of grey system models such as GM (1, 1) model with DGM (1, 1) [42].

Conclusions
This study employed the proposed W-AL, single ARIMA and single LSTM to predict precipitation, using monthly data over a long period of from1967 to 2017 for three stations. Additionally, drought analysis by using the CZI drought index, using annual precipitation amounts for the same years was carried out. Finally, drought years, as classified by the CZI, applying the grey prediction models of GM (1, 1) and DGM (1, 1), were predicted. The following main conclusions are drawn from this study.
The proposed W-AL model at different ratios of training and test sets all exhibited higher prediction accuracy than the ARIMA and LSTM, based on different climate types for monthly precipitation data. In addition, by comparing the R2 values obtained by the W-AL models of the three stations, it can be found that the fitting effect of the W-AL method in Linjiang with humid and semi-humid climate type was best and was followed by Changchun with semiarid and semihumid climate type and Qian Gorlos with arid and semiarid type.
The drought index CZI results revealed that, drought events have become more frequent since 1987, and that all stations experienced some levels of drought in the 1980s, mid-2000s and early 2010s. On the other hand, the results indicated that there were more numerous slight drought events in Changchun and more numerous moderate drought events in Linjiang and Qian Gorlos.
GM (1, 1) and DGM (1, 1) were used to predict drought years. According to the results, GM (1, 1) always showed higher accuracy than DGM (1, 1) at different climate types, with an average relative error of 2.22% at a minimum and 6.66% at a maximum. Therefore, GM (1, 1) was applied to predict drought years that were close to the actual conditions. Additionally, the best prediction effect for drought events was relatively arid areas, and the worst was a relatively humid area.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
Publicly available datasets were analyzed in this study. This data can be found here: http://data.cma.cn/.