Volatility Analysis and Visualization of Climate Data Based on Wavelets

In this article, we analyze the real meteorological data recorded by Wenzhou Meteorological Bureau from 1951 to 1997. The data has not been used elsewhere and is available at Meteorological Station Wenzhou (ID: CHM00058659) at https://geographic.org/global weather/china. We perform the time series volatility analysis including ARMA, ARIMA, ARCH-LM, PARCH, SARMA and Morlet wavelet analysis and use the Mann-Kendall (M-K) test to analyze both the trend and mutation defined by statistics sequence. In addition, a Morete wavelet time-frequency model is established to show that both the precipitation and temperature have a very important 12-month cycle and the precipitation is also very unstable. We then employ the STL, coif1 decompositions and NAR model to capture both the volatility and Heteroscedasticity in the data. In addition, the performance of the fitted model has been proven to be satisfactory on actual climate data with the small Mean Square Error (MSE), Root-Mean-Squarred Error (RMSE), and coefficient of determination. Finally, monthly average temperature is added as an exogenous (covariate) variable and a nonlinear autoregressive exogenous model is employed to improve the performance of the model. Our results show that the performance of NARX model is more accurate and stable with better mean square error.


Introduction
Climate change is an important environmental, social and economic issue. It threatens the achievement of Millennium development goals aimed at poverty and hunger reduction, health improvement and environmental sustainability. A large number of studies have been conducted on several aspects to address the impact of climate change. Rosenzweig and Parry (Rosenzweig & Pary, 1994) provides a global assessment of climate change on food supply and predict grain yield losses in several countries. (Barrios, Quattara & Strolbl, 2008) recommend using the regression analysis to estimate the impact of weather changes on agriculture in some sub-Saharan African countries.
One of the important factors to measure climate change is to observe changes in temperature and precipitation in a specific area. In the past few decades, different research groups have released different global and hemispheric temperature estimates based on existing observational data. Significant improvements in information recovery and processing analysis techniques have led to longer temperature records based on more accurate estimates and better uncertainty assessments (Hansen et. al., 2010). Due to their importance, many forecast models for weather forecasting have been developed. Time series analysis has a wide range of applications in engineering, economics, meteorology, finance and many other fields. A number of stochastic time series models such as the Markov, Box-Jenkins (BJ) Seasonal Autoregressive Integrated Moving Average (SARIMA), deseaonalized Autoregressive Moving Average (DARMA), Periodic Autoregressive (PAR), Transfer Function Noise (TFN) and Periodic Transfer Function Noise (PTFN), are in use of these purposes (Ediger & Akar, 2007). In recent years, a review work related to the detection and attribution of climate change based on time series techniques is presented (Estrada & Perron, 2014) and the authors used co-integration techniques to solve the problem and addressed the point that time series can be better described as the stationary fluctuations of the trend function with the slope. In addition, some univariate time series models that deal with only one time series, including the autoregressive integrated moving average (ARIMA) model and its derivatives such as seasonal ARIMA (SARIMA), periodic ARIMA, and deseasonalized ARIMA models, have long been applied in monthly streamflow forecasting (Bender & Simonovic, 1994).
The features of variable and their interpretations are illustrated in Table 1. Note that there are total of 564 observations recorded and some missing data are found in PRCP, TAVG, TMAX amd TMIN. Given that the ratio of missing data to the total observations is small, so we choose to utilize the Catmull-Rom Spline method to smooth and fit the missing data.
The summary of basic statistics as well as numerical attributes of each variable are shown in Table 2, where the mean, median, min, and max are self-explanatory variables. Figs. 1-2 show the time series of featuring variables, namely PRCP, TAVG, TMAX and TMIN, in the scale of Month. Here one can see that the fluctuations of PRCP are less cyclical than that of other variables. However, the volatility is stronger compared with other variables, indicating the nature of time series of precipitation. Notice that the trend of time series of TAVG, Tmax and Tmin are basically similar, and all of the characteristic variables appeared in obvious annual cycles. In other words, the PRCP fluctuates around 5, while TAVG, TMAX and TMIN are relatively maintained at a certain temperature level of 65 • F in the scale of year.
In Fig. 3, it is seen that the PRCP started to increase slowly from January, and data numbers reach the peak both in June and August and decrease in the other months, showing a shape of zigzag. It can also be observed that TAVG, TMIN and TMAX have a similar trend in the time series, and the overall change is slow. In addition, the temperature began to rise slowly in January, peaked in July and August, and then decreased slowly in the next few months until it reaches the lowest point in January and February.
In order to make the comparisons of featuring variables of TAVG and PRCP in time series analysis, we select the top 10 years of values in July and last 10 years of values in January from TAVG, while for PRCP, we select the June data of the past 10 years and the January data of the last 10 years to facilitate the model processing. The results are shown in Fig. 4. Fig. 4 shows that the highest temperature in July occurred in 1988, and the lowest temperature in January occurred in 1963. In addition, it also shows that the maximum precipitation occurred in June of 1990 and the minimum precipitation in occurred in January of 1963. It is noted that both the minimum temperature and precipitation occurred in January of 1963, and the ten highest temperatures in July in history include 1963, pointing out that extreme weather occurred in 1963.
The scatter plots in Fig. 5 shows that TAVG, TMIN, TMAX are all highly linear correlated, which confirmed the results obtained in the Figs. 1-2. Here Fig. 5 also gives us the idea to choose TAVG variable to represent the temperature to avoid the Pseudo-regression.

Periodicity Analysis
Periodicity occurs naturally in many environmental time series, such as hourly tide levels, daily stream flows, monthly average temperatures, and carbon dioxide exchange by growing plants. This section will focus on the examples of monthly precipitation and monthly average temperature in a demonstration of methods for analyzing periodically correlated time series together with the use of the Mann-Kendall (M-K) test and Morelet wavelet analysis.

Mann-Kendall Test
It is seen from previous sections that the time series of year trend of PRCP is fluctuating, therefore it is essential to rule out the Mutations of the series. In order to examine the mutational points, the Mann-Kendall (M-K) test is carried out on the time series data. The M-K method is a nonparametric statistical test method and is also known as a No-Distribution test. This method does not need to set the data distribution characteristics in advance, and is not interfered by some missing values and outliers. Therefore it is suitable for trend test and mutation points of meteorological data and hydrological time series (Yue, Pilon & Cavadias, 2002). In addition, M-K test can be used to analyze both the trend and mutation of the time series by defining statistics values of UF and UB. Here the UF is defined as where and On the other hand, the UB is defined as By analyzing the statistical sequence of UF and UB, the trend and change of time series can be analyzed, and the time and area of mutation points can be identified. If UB > 0, the time series is indicating an upward trend, while if UF < 0, the time series indicates a downward trend over time. When it comes to exceeding the critical line, the data numbers indicate either a significant upward or downward trend. If both the time series curves of UF and UB have intersection points and the intersection point is between the critical lines, then the intersection point is the time when the mutation begins. The results of testing the time series curves of UF and UB are shown in Fig. 6. It can be seen in Fig. 6 that both the statistic values of UF and UB are developed within two critical line, which means there are no mutations between 1951-1997.

Monthly/Yearly Periodicity Analysis Based on Morlet Wavelet
It can be seen from the Figs. 1-2 that both the PRCP and TAVG have significantly Monthly period characteristics from 1951 to 1997. However, in the "year" scale, it can be seen that the period characteristics are not as important as the "month" scale, especially during this period, the TAVG remains almost unchanged. In order to quantitatively detect the periodic nature of data, the wavelet time-frequency (Morlet wavelet) method is employed to process and analyze the time series data. The Morlet wavelet is a continuous plane wave (Torrence & Compo, 1998) modulated by Gauss function based on the following function: The formula of discrete wavelet transform is expressed as: where * denotes the complex conjugate, a is a scale factor related to period and frequency, b is translation factor of time position, i is the time position label of data series, and f (t) is the variable of time series. In addition, W f (a, b) is a wavelet coefficient and δt is the difference of time series. The wavelet power spectrum is defined as follows: The total wavelet power spectrum E a represents the energy density corresponding to different scales a, which is defined as In addition, it is found that the significant test of wavelet power spectrum is quite vital, which could be used to test the standard spectrum of red noise or white noise. The test process follows the patter: calculate the effective degree of freedom of wavelet power spectrum distribution, carry out the significance test of χ 2 distribution. As such, the theoretical power spectrum P of red noise or white noise will be calculated accordingly. If the total wavelet power spectrum (E a ) is larger than the theoretical spectrum, the test is shown to be significantly credible. Here the theoretical power spectrum is defined as: where χ 2 ν is the value of χ 2 with degree of freedom ν, and σ 2 is the variance of time series which is defined by Then the theoretical spectrum could be calculated by the following formula: where α is the the Autocorrelation coefficient of lag 1 of series. if E a > P a , then it shows that the period corresponding to the total wavelet power spectrum is significant. To reduce side effects, the original time series was extended by 10 values from both sides.
Figs. 7-8 shows the monthly average precipitation in Wenzhou, eastern China from 1951 to 1997, which contains several periodic variations in different scales, forming oscillation centers between positive and negative scales with obvious interannual variations. It is observed that there are mainly three distinct periodic changing times in the process, namely the 12 months, 30-50 months and 80-120 months. However, from the analysis of wavelet power spectrum (Fig. 7, right color bar), it is found that only the 12-month period time change passes the noise test (see Fig. 7, red dotted line). From the time series of 12-month (one-year) analysis, it is found that the monthly precipitation periodic oscillation in Wenzhou area is very significant, and there are obvious clustering oscillations appeared in the whole time domain. It can also be seen that there are five clustering oscillations, namely 1951-1955; 1959-1963, 1968-1972, 1976-1980, and 1985-1989, appeared in the model which can be confirmed in Fig. 8. The clustering oscillation is manifested as the clustering oscillation of Real Part. In addition, we can know from the wavelet power spectrum that the wavelet power values in the early (1951)(1952)(1953)(1954)(1955)(1956)(1957)(1958)(1959)(1960)(1961) and late (1986)(1987)(1988)(1989)(1990)(1991) periods are larger than those in the middle . Notice that there are two peaks in the year of 1986 and 1991, which indicate that the 12-month cycle changing time of precipitation in these two specific years is very significant.
Figs. 9-10 show the average monthly temperature in Wenzhou from 1951 to 1997, which contains almost only one scale of periodic change (about 12 months) and notice that its wavelet power spectrum passed the noise test (red dotted line is red noise spectrum).The annual cycle oscillation of monthly temperatures in Wenzhou area is very significant throughout the whole time domain. The results in Fig. 10 show that there are similar periodic changes in the entire time domain, which shows that the annual periodic change of monthly average temperature is more significant than that of monthly average precipitation, and the interference is small.

Volatility Analysis On PRCP
The understanding of volatility is a key concept for analyzing environmental data. This is of greater importance for climate data because it provides key aspects such as the return on precipitation and temperature of weather conditions and helps with effective analysis. The unpredictable nature of volatility causes heteroskedasticity which leads to difficulty in modelling. Consequently, time series models are desirable to predict volatility. An illustration of the same has been shown through an example of fitting time series models on the volatility of a listing from the the Stock Exchange. This paper attempts to employ different methods to analyze the volatility phenomenon. Based on the results shown in Fig. 8, it is seen that the spectrum is not significantly zero besides the sole peak. It is also observed that there is a flat peak between 80-100 months and the spectrum remains at a non-zero level, which means that although PRCP fails to pass the red noise test, it still has a certain degree of volatility. Those volatility might play a vital role in dominating the regularity or uncertainty of precipitation in wenzhou. Therefore, it is necessary to analyze it in this article. The flow chart of research process is demonstrated in Fig. 11 The remainder of the time series is extracted based on Seasonal and Trend decomposition using Loess (STL) procedure. In order to capture the volatility in remainder analysis, the technique of wavelet decomposition is employed to decompose the time series into low-frequency and high-frequency terms of d1, d2, d3, d4 and d5 through the method of wavelet decomposition and reconstruction. The time series models ARMA, SARMA-PARCH, SARMA-EGARCH, and NAR are established according to their respective characteristics. Finally, we fit the values linearly and combine the values to obtain the final fitting of R series.

STL Decomposition
STL stands for Seasonal and Trend decomposition using Loess. This is a statistical method of decomposing a Time series data into three components containing seasonality, trend and residual, whereas seasonality is a regular and predictable pattern that recur at a fixed interval of time.
In order to analyse the volatility of PRCP, some certain factors should be ruled out. Given that there is 12-months periodic phenomenon in PRCP, we use STL procedure to process data. Based on STL procedure (Robert, William & Irma, 1990), PRCP data will be decomposed into three parts, namely seasonal, trend and remainder. The PRCP data can be represented as follows: where S t and T t are terms meaning certainty, and the term R t is the remainder due to the fact that there might contain some vital information about the volatility. Based on the previous analysis and monthly data of precipitation and temperature, we set the periodicity to 12 and smooth the windows of seasons, trends and filters. These are the default functions in the software EVIEWS, and the results in seasonal and diagnostics panels are shown in Figs. 12-13, suggesting that seasonal component is around 12 months and the trend component has no significantly upward or downward trend. Flat lines in each of months indicate that the choice of 12 months as the periodicity and parameters for seasonal window is appropriate.

Coif1 Wavelet Decomposition
It is well known that the wavelet analysis method has a strong multi-resolution ability. It can separate the high frequency (representing the period and disturbance part) and the low frequency (representing the trend part) of the precipitation sequence, and process the PRCP sequence from different levels according to the characteristics of each sub-sequence. In this section, we resort to the coif1 wavelet analysis. The results of decomposition together with residuals of reconstruction are shown in Figs. 12-13. In addition, Figs. 14-15 illustrate the coif1 wavelet decomposition of PRCP data and Fig. 15 shows that the residual is about 10 −11 order, which indicates the reconstruction is efficient.

Models
In this section, different models are selected for variables a5, d1, d2, d5 from Table 7-10. First,the ARMA model requires the property of sequence stationarity. In general, time series data are prone to non-stationary problems which can lead to the destruction of the consistency of large sample statistics and produce pseudo regression problem. Before the follow-up test, the stationarity of each variable is tested by the augmented Dickey-Fuller test (Dickey & Fuller, 1979), which is used to test the null hypothesis that a unit root is present in an autoregressive model. The test results are shown in Table 3.
Here D (*) denotes the first-order difference, and [C, T, L] is the basic type of unit root, which represents constant term, trend term and lag order term respectively. The probabilities in [**] represent one-side p-values suggested by Mackinnon (1996). In addition, results shows that six variables are stationary at 5% level. The lags of ARMA /SARMA model are selected based on minimum AIC. The procedure of model selection and its criteria based on AIC and dependent variables a5, d1, d2, and d5 can be seen in Table 7-10. It is the selection of parameters of ARMA(p,q) based on AIC, SC and HQ criteria. In addition, some variance equations are added in order to capture the volatility and conditional variance. If the residual has shown to exist in Heteroscedasticity, it means that it is possible to establish models to fit the data in the sense of Heteroscedasticity. In general, heteroscedasticity of residuals could be tested based on ARCH-LM model by setting a optimal lags operator, and the choice of the best lag factor varies by methods. In (Chong, Ahmad & Abdullah, 1999), the authors use ARCH-LM(12) to reject the null hypothesis and concluded that a very high-order ARCH model is needed to model the heteroscedasticity. In (Engle, 1982), it is shown that the Lagrange multiplier test for a first-order linear ARCH effect is not significant. However, testing for a fourth-order linear ARCH process and the chi-squared statistic with 4 degrees of freedom will be highly significant. In (Engle, 2001), authors choose one lag in order to incorporate the ARCH effect. Given that the ARCH-LM testing might be sensitive to the choice of lags which means different lags could give different results and the auto-correlation in most series dies out after 15 lags, we choose the number of lags 1,5,10,15 for research, and choose the best lag according to the AIC, SC and HQ criteria. We choose the lag with the smallest number of AIC, SC, HQ to judge heteroscedasticity.
In addition, the ARCH model family is considered an effective way to solve our problems. There are lots of research results addressing the heteroscendasticity problem of the data. In (Zhou & Li, 2018), the authors use the Generalized Auto-Regressive Conditional Heteroskedasticity (GARCH) model to analyze the carbon price fluctuation characteristics.
In (Amiri, 2014), the authors analyze the conditional variance of a function of ACDC level known as ACDC level growth rate (ACDCGR) using the generalised autoregressive conditional heteroskedasticity (GARCH) and GARCH models with leverage effect. In addition, research results in (Liu & Shi, 2013) show that SARMA-GARCH,ARMAGARCH, ARMA-EGARCH-M models can be used to deal with different volatility problems, which shows the efficiency and effectiveness of ARCH model family. In this article, through the time-testing and selections based on the criteria information, residual diagnostic and stability of model, we utilize the PARCH and EGARCH as variance equations to capture the volatility for variables d1, d2, and d5. The variance equation of PARCH(p,q) is expressed as follows: where δ > 0, |γ i | ≤ 1 for i = 1, . . . , r, with γ i = 0 for all i > r, and r ≤ p. Moreover, the variance equation of EGARCH(p,q) is: Notice that the EGARCH model first relaxed the non negative constraints on the parameters of the GARCH model, and then adds an asymmetric term to represent the leverage effect. In our approach, the ARMA model is employed for a5 due to the fact that there is no heteroskedasticity in the residual analysis. However, as for d3 and d4, we find that there are no suitable SARMA and ARCH models for our data analysis. The main problem is that most models can not pass the residual white noise test. It is believed that the model parameters d3, d4 could be highly complex and non-linear. Therefore, nonlinear autoregressive (NAR) based on neural networks will be one of our methods to process the data in the following sections. In this study, the calculated mean square error MS E and the coefficient of determination R 2 are used to compare the performance of the model.
where n is the number of data points, y i is observed values,ŷ i is predicted values, andȳ is mean of data y.
The processes of selecting the parameters of d1, d2, d5, and a5 from Tables 7-10 are similar. Here we demonstrate the selection of d1 as an example. Firstly, the mean equation is selected based on AIC (see Tables 7-10), and from that we choose ARMA(1,2) to process the data. In what follows, the ARCH-LM test is employed for testing the Heteroscedasticity. The results are shown in Table 4, which shows that the null hypothesis is rejected in all lags. Based on AIC, SC, HQ, we choose lag10 from the table to be the optimal lag with P value set to 0, indicating that there is Heteroscedasticity in residuals ( d1 t ). Then the squared-residuals-autocorrelation is tested and results are shown in Fig. 16. The results suggest that the squared residuals shows autocorrelation of up to lag 15, indicating that it is not white noise and a suitable model can be established to extract the information. It also indicates that Heteroscedasticity model could be built using the residuals ( d1 t ). Through times of testing and selection based on criteria information, residual diagnostic and stability of the model, we finally choose ARMA(1,2)-PARCH(2,1) as a good candidate model to perform the data. Its equation is expressed as ( where L is the lag operator and the upper index d1 stands for the variable name . The parameter estimations for variables d1, d2, d3, d4 and d5 shown in Tables 11-14 demonstrate that the parameters in mean equations are all significant, while the parameters in variance equation are significant apart from γ 1 . It also indicates that there is no significant Leverage effect. Table 5 shows the results of ARCH-LM test after building the PARCH model. Based on the optimal information criteria (AIC, SC and HQ), lag1 is selected by the model and p value is 0.8498, which indicates that there is no Heteroskedasticity in the residual. The residuals autocorrelation diagrams in Figs. 23-28 suggests that there is no residuals autocorrelation at 5% level. It is suggesting that the model is effective and can be used to fit the data well. The fitting results and conditional variance diagram are shown in Fig. 17 where the red dot line is standard error based on variance equation. The variance diagram shows that uncertainty volatility of d1 is the greatest in data source No.1974-M01 and smallest in data source No.1961M07-1969M12.
In a similar fashion, the SARMA(5,5)(1,0)-EARCH(1,1) models can be established for d2 with the following resulting ( In a similar fashion, the SARMA(3,4)(1,1)-PARCH(1,1) model can be established for d5 with the resulting equations as: To proceed, the ARMA(4,4) model can be established for a5, and the resulting equations are: Notice that the NAR model can be established for d3, d4 in a similar fashion. The model could be implemented by neural network. Therefore we create a two layers with 15 hidden neurons and set the numbers of delay to be 12 based on previous periodicity analysis.  . 20) and the residuals autocorrelation graphs in Figs. 23-28 also show the residuals are white noise and the information has been fully extracted. Notice that the total fitted remainderR t of the model can be represented as linear combination of all sub-series as follows:R t =â 5 t +d 1 t +d 2 t +d 3 t +d 4 t +d 5 t .
Here we set the fitted seriesR t from the real Data Source: No.1955M04-1997M12 due to the fact that the lags in some ARMA models may cause some missing data before 1955M04. Consequently, the result shows that MSE is 2.5008 and RMSE is 1.5814,R 2 =0.8057, which indicates that 80% of the changes in volatility can be interpreted by model.

Covariates and NARX Model
In order to compare and improve the model that can accurately fit the volatility, we consider the feature variable of TAVG as exogenous variable and put it into our models for test. In view of the obvious periodicity of TAVG, it cannot be regarded as a strictly stationary sequence. Therefore, it is not suitable to use the ARIMA model. In general, this is not a fixed sequence and the use of ARIMA model requires that each sequence is stable, otherwise there will be a problem of estimation in parameters. As a result, we proposed using a nonlinear autoregressive exogenous (NARX) neural network model to convert the input into two neurons so that TAVG can be added in the model. The NARX can be described by the following non-linear function: where x(t−d) is an exogenous variable, and d is the time delays which can be used to determine the effectiveness of NARX. In time series modeling, a nonlinear autoregressive exogenous model (NARX) is a nonlinear autoregressive model which has exogenous inputs. This means that the model relates the current value of a time series to both the past values of the same series and current and past values of the driving (exogenous) series-that is, of the externally determined series that influences the series of interest. The NARX model has a wide range of applications in multi-step neural network prediction and shows a good performance in dealing with some complex series that can not be fitted well by ARIMA.
In (Abou Rjeily et. al., 2017), NARX was built to develop a flooding forecast system and showed a high efficiency in representing the real urban drainage systems response to a forecasted storm event in flooding forecast system. In (Guzman, Paz & Tagert, 2017), the authors employed NARX to simulate daily groundwater levels at a local scale in the Mississippi River Valley Alluvial aquifer located in the southeastern United States. In (Wunsch, Liesch & Broda, 2018), NARX was applied to obtain groundwater level forecasts for several wells in southwest Germany and showed an outstanding suitability in model prediction for groundwater level. In this section, we construct a NARX Neural network as shown in Fig. 21 and use the Bayesian Regularization algorithm The NARX network to train the proposed model. As a result, the optimal parameters including the number of hidden neurons and time delays are selected based on the performance of Residuals auto-correlations MS E and R 2 . Firstly, each variables such asâ 5 t ,d 1 t ,d 2 t ,d 3 t ,d 4 t ,d 5 t are fitted by NARX model, and then calculatesR t by applying the equation in (21). The results are illustrated in Table 6. It shows that the MSE is 0.5219, which is smaller than the model case without exogenous variables (MSE=2.5008). In addition, it also shows that R 2 is 0.9195, accounting for nearly 92% of the change in volatility explained by the model. The result of fitted data from a new model is shown in Fig. 22. In our view point, the result of modified model with R 2 = 0.9195 is better than previous model without exogenous variable (R 2 = 0.8057).