# Forecasting Uranium Resource Price Prediction by Extreme Learning Machine with Empirical Mode Decomposition and Phase Space Reconstruction

Authors: Qisheng Yan; Shitong Wang; Bingqing Li
Journal: Discrete Dynamics in Nature and Society
 (2014)

---

## Abstract

A hybrid forecasting approach combining empirical mode decomposition (EMD), phase space reconstruction (PSR), and extreme learning machine (ELM) for international uranium resource prices is proposed. In the first stage, the original uranium resource price series are first decomposed into a finite number of independent intrinsic mode functions (IMFs), with different frequencies. In the second stage, the IMFs are composed into three subseries based on the fine-to-coarse reconstruction rule. In the third stage, based on phase space reconstruction, different ELM models are used to model and forecast the three subseries, respectively, according to the intrinsic characteristic time scales. Finally, in the foruth stage, these forecasting results are combined to output the ultimate forecasting result. Experimental results from real uranium resource price data demonstrate that the proposed hybrid forecasting method outperforms RBF neural network (RBFNN) and single ELM in terms of RMSE, MAE, and DS.


---

## Body

## 1. Introduction

Uranium resource products have been widely used in economics, military, social lives, and so on and have a revolutionary effect on the real world in many areas. Uranium resource is both the material basis for the development of nuclear energy and a kind of strategic resource. Therefore, more accurate forecasts for international uranium resource prices play an increasingly important role in the development planning and utilization of nuclear energy.Forecasting the uranium resource prices is one of the most important and challenging tasks due to its inherent nonlinearity and nonstationary characteristics. In the past decades, price prediction has attracted increasing attention by lots of academic researchers. The forecasting approaches used in the literature can be classified into two categories: statistical models and artificial intelligence models [1, 2]. However, the statistical models cannot effectively capture nonlinear patterns hidden in price time series owing to the fact that these models are developed based on the underlying assumption that the time series being forecasted are linear and stationary [3]. In order to overcome this limitation of statistical models, a good deal of nonlinear models have been proposed, among which the artificial neural network (ANN) has attracted a growing interest by researchers due to their excellent nonlinear modeling capability [4–8]. Many studies conclude that the ANN model outperforms various conventional statistical models. However, ANN suffers from local minimum traps and the difficulty of determining the hidden layer size and learning rate [9]. A new learning algorithm for the single-hidden-layer feedforward neural network (SLFN) called the extreme learning machine (ELM) has been proposed recently and overcomes the aforementioned disadvantages [10, 11]. In the learning process of ELM, the input weights and hidden biases are randomly chosen, and the output weights are analytically determined by using the Moore-Penrose generalized inverse. ELM can learn much faster with a higher generalization performance than the traditional gradient-based learning algorithms and solves the problem of stopping criteria, learning rate, learning epochs, and local minima [10–13]. In recent years, ELM has attracted a lot of attention and has become an important method in nonlinear modeling [11–13].When building intelligent prediction models directly using original values, it is difficult to obtain satisfactory forecast results due to the high-frequency, nonstationary, and chaotic properties of uranium resource price data. Hence, in order to further improve the prediction performance, before constructing a forecasting model, recent research efforts on modeling time series with complex nonlinearity, dynamic variation, and high irregularity are to initially utilize information extraction techniques to extract features hidden in the data, then use these extracted characteristics to construct the forecasting model [14–18]. That is to say, by some means of suitable feature extractions or signal processing methods, the useful or interesting information which may not be observed directly from the original data can be revealed in the extracted features. Thereby, an effective forecasting model possessing better prediction precision will be developed.Empirical mode decomposition (EMD), based on Hilbert-Huang Transform (HHT), is very suitable for decomposing nonlinear and nonstationary time series, which adaptively represents the local characteristic of the given signal [19, 20]. By using EMD, any complicated signal can be decomposed into a finite and often small number of Intrinsic Mode Functions (IMFs), which have simpler frequency components and stronger correlations, thus are easier and more accurate to forecast [8]. Recently, the EMD has been widely used in many fields, such as in the analysis of the atmosphere time series [21], river water turbidity forecasting [22], crude oil price prediction [23], short-term wind power prediction, and so forth, [8, 13, 14, 16, 24].In this study, we propose a hybrid uranium resource prices forecasting model by integrating EMD, phase space reconstruction (PSR), and extreme learning machine (ELM). Firstly, the original uranium resource price series are first decomposed into a finite number of independent intrinsic mode functions (IMFs), with different frequencies. Secondly, the IMFs are composed into three subseries based on the fine-to-coarse reconstruction rule. And then, based on phase space reconstruction, different ELM models are used to model and forecast the three subseries, respectively, according to the intrinsic characteristic time scales. Finally, these forecasting results are combined with the ultimate forecasting result output. Moreover, experimental results from real uranium resource price data demonstrate that the proposed hybrid forecasting method outperforms RBF neural network (RBFNN) and single ELM in terms of RMSE, MAE, and DS.The rest of this paper is organized as follows. Section2 gives brief overviews of EMD, PSR, and ELM. The proposed model is described in Section 3. Section 4 compares the experimental results obtained by the proposed hybrid approach and those by RBF neural network and single ELM, and this paper is concluded in Section 5.


## 2. Methodology

### 2.1. Empirical Mode Decomposition (EMD)

EMD is a new signal processing technique. Unlike wavelet decomposition, EMD is not required to determine a filter base function before decomposition. The main idea of EMD is to decompose original time series data into a sum of oscillatory functions, namely, intrinsic mode functions (IMFs). In the EMD, the IMFs must satisfy two conditions:(1) the number of extrema (sum of maxima and minima) and the number of zero crossings must either equal or differ at most by one, and (2) at any point, the mean value of the envelope defined by the local maxima and the envelope defined by the local minima is zero.The essence of EMD is the sifting process which extracts IMFs from the original data. The algorithm of EMD is described as follows [19, 20, 25].Step 1.
Identify all the local extrema including the minimum values and maximum values in time series datax(t).Step 2.
Generate the upper and lower envelopesemax⁡(t) and emin⁡(t) by a cubic spline line.Step 3.
Calculate the mean valuem1(t) from the upper and lower envelopes and then generate the mean envelope as
(1)m1(t)=[emax⁡(t)+emin⁡(t)]2.Step 4.
Calculate the difference between the time series datax(t) and the mean value m1(t). The first difference h1(t) is designed as a protointrinsic mode function
(2)h1(t)=x(t)-m1(t).Step 5.
Check whether the protointrinsic mode functionh1(t) satisfies the properties of IMF or not. Ideally, h1(t) should be an IMF. However, it may generate a new extremum and shift or exaggerate the existing extrema in the sifting process.
If properties ofh1(t) satisfy all the requirements of an IMF, h1(t) is denoted as the ith IMF ci(t) and substitutes the residue r1(t) for the original time series data x(t); that is,
(3)r1(t)=x(t)-h1(t).
Otherwise, h1(t) is not an IMF. Then, it substitutes h1(t) for the original time series x(t).Step 6.
Repeat from Step1 to Step 5. The sifting process stops when the residue satisfies one of the termination criteria. First, the residue or the ith component is smaller than the predetermined threshold or becomes a monotonic function such that no more IMF can be extracted. Second, the number of zero crossings and extrema is the same as that of the successive sifting step.By using the above algorithm, the original time series datax(t) can be decomposed into n modes and a residue as follows:
(4)x(t)=∑i=1nci(t)+rn(t),
where n is the number of IMFs, ci(t) represents IMFs which are nearly orthogonal to each other and periodic, and rn(t) is the final residue which is a constant or a trend. By the sifting process, each IMF is independent and specific for expressing the local characteristics of the original time series data. The set of IMFs is derived from high frequency to low frequency, while rn(t) represents the central tendency of data series x(t). In addition, EMD can also be taken as a filter of high pass, band pass, or low pass.


### 2.2. Phase Space Reconstruction (PSR)

Takens embedding theorem [26] provides theoretical foundation for the analysis of time series generated by nonlinear dynamical systems. Later, Sauer et al. [27] show that a phase space can be reconstructed from a univariate chaotic time series. Let a univariate time series {xk}k=1n, where N is the length of the time series, generated from a d-dimension chaotic attractor, then a phase space Rd of the attractor can be reconstructed by using the delay coordinate defined as
(5)x(k)=[xk-τ,…,xk-mτ],
where m is known as the embedding dimension of reconstructed phase space and τ is the delay constant.The selection of the embedding dimensionm and the delay constant τ is very important for prediction modeling [3]. Therefore, for a given delay time τ, a time series {xk}k=1n is represented in the so-called “phase space” by a set of delay vectors (DVs) x(k)=[xk-τ,…,xk-mτ] of a given embedding dimension m.(
1) Determination of the Delay Time
τ
and the Minimum Embedding Dimension
m. The entropy ratio (ER) method [28] is a novel method for determining the set of parameters for a phase space representation of a time series. Based upon the differential entropy, both the optimal embedding dimension and time lag are simultaneously determined.Based upon the probability density function of data, the differential entropy is defined as(6)H(x)=-∫-∞+∞p(x)ln⁡p(x)dx.Particularly convenient is the Kozachenko-Leonenko (K-L) estimate of the differential entropy(7)H(x)=∑j=1Nln⁡(Nρj)+ln⁡2+CE,
where N is the number of samples in the dataset, ρj is the Euclidean distance of the jth delay vector to its nearest neighbour, and CE is the Euler constant. For a given embedding dimension m and time lag τ, let H(x,m,τ) denote the differential entropies estimated for time delay embedded versions of a time series x.The K-L estimates for the time delay embedded versions of the original time seriesH(x,m,τ) and H(xs,i,m,τ) are computed using (7) for increasing m and τ (index i refers to the ith surrogate). To determine the optimal embedding parameters, the ratio
(8)I(m,τ)=H(x,m,τ)〈H(xs,i,m,τ)〉i
needs to be minimized, where 〈·〉 denotes the average over i. To penalize for higher embedding dimensions, the minimum description length (MDL) method is superimposed, yielding the “entropy ratio” (ER)
(9)Rent(m,τ)=I(m,τ)+mln⁡NN,
where N is the number of delay vectors, which is kept constant for all values of m and τ under consideration. The minimum of the plot of the entropy ratio yields the optimal set of embedding parameters.(
2) Identification of Nonlinearity. The delay vector variance (DVV) method [29] is a novel analysis of a time series which examines a signal’s unpredictability by observing the variability of the targets belonging to sets of similar delay vectors (DVs). The DVV method can be summarized as follows for a given embedding dimension.(i)
The meanμd and standard deviation σd are computed over all pairwise distances between DVs.(ii)
The setsΩk are generated, which consist of all DVs that lie closer to x(k) than a certain distance. The distances are taken from the interval [μd-ndσd,μd+ndσd], for example, uniformly spaced, where nd is a parameter controlling the span over which to perform the DVV analysis.(iii)
For every setΩk, the variance of the corresponding targets, σk2, is computed. The average over all sets, divided by the variance of the time series σx2, yields the measure of unpredictability σ*2:
(10)σ*2=(1/N)∑k=1Nσk2σx2.In the following study, the linear or nonlinear nature of the time series is examined by performing DVV analyses on both the original and a number of surrogate time series, using the optimal embedding dimension of the original time series.(
3) Identification of Chaotic Characteristic. When the largest Lyapunov exponent of the system is larger than zero, it indicates that there is a chaotic attractor which can be used to measure the chaotic degree [30]. The largest Lyapunov exponent is computed by the Wolf method [31].


### 2.3. Extreme Learning Machine (ELM)

Suppose there areN distinct samples {(xk,tk)}k=1N, where xk∈Rn and tk∈Rm. The SLFN with M hidden neurons can be described as
(11)∑i=1Mβig(wi·xk+bi)=ok,k=1,2,…,N,
where wi=[ωi1,ωi2,…,ωin]T  is the weight vector connecting the ith hidden neuron and the input neurons, βi=[βi1,βi2,…,βim]T  is the weight vector connecting the ith hidden neuron and the output neurons, oi=[oi1,oi2,…,oim]T  is the actual output vector, and bi is the threshold of the ith hidden neuron. g(·) represents the activation function of hidden neuron, and the operation wi·xk denotes the inner product of wi and xk.If the SLFN can approximate theN samples with a zero error then we have
(12)∑k=1N∥ok-tk∥=0.Thus, there also exist parametersβi, wi, and bi such that
(13)∑i=1Mβig(wi·xk+bi)=tk,k=1,2,…,N.The above equations can be compactly described as(14)Hβ=T,
where
(15)H(w1,w2,…,wM,b1,b2,…,bM,x1,x2,…,xN)=[g(w1·x1+b1)⋯g(wM·x1+bM)⋮⋮⋮g(w1·xN+b1)⋯g(wM·xN+bM)]N×M,β=[β1T⋮βMT]M×m,T=[t1T⋮tNT]N×m.Unlike the traditional function approximation theories which require the adjustment of input weights and hidden layer biases, the input weights and hidden biases are randomly generated. Thus, training an SLFN is simply equivalent to finding a least-squares solution of the linear functionHβ=T. Consider
(16)∥Hβ^-T∥=min⁡β∥Hβ-T∥.
The smallest norm least-squares solution of the above linear system is
(17)β^=H+T,
where H+ is the Moore-Penrose generalized inverse of matrix H. Owing to the Moore-Penrose generalized inverse, the learning speed is dramatically increased for the single hidden-layer feedforward neural network [12].


## 2.1. Empirical Mode Decomposition (EMD)

EMD is a new signal processing technique. Unlike wavelet decomposition, EMD is not required to determine a filter base function before decomposition. The main idea of EMD is to decompose original time series data into a sum of oscillatory functions, namely, intrinsic mode functions (IMFs). In the EMD, the IMFs must satisfy two conditions:(1) the number of extrema (sum of maxima and minima) and the number of zero crossings must either equal or differ at most by one, and (2) at any point, the mean value of the envelope defined by the local maxima and the envelope defined by the local minima is zero.The essence of EMD is the sifting process which extracts IMFs from the original data. The algorithm of EMD is described as follows [19, 20, 25].Step 1.
Identify all the local extrema including the minimum values and maximum values in time series datax(t).Step 2.
Generate the upper and lower envelopesemax⁡(t) and emin⁡(t) by a cubic spline line.Step 3.
Calculate the mean valuem1(t) from the upper and lower envelopes and then generate the mean envelope as
(1)m1(t)=[emax⁡(t)+emin⁡(t)]2.Step 4.
Calculate the difference between the time series datax(t) and the mean value m1(t). The first difference h1(t) is designed as a protointrinsic mode function
(2)h1(t)=x(t)-m1(t).Step 5.
Check whether the protointrinsic mode functionh1(t) satisfies the properties of IMF or not. Ideally, h1(t) should be an IMF. However, it may generate a new extremum and shift or exaggerate the existing extrema in the sifting process.
If properties ofh1(t) satisfy all the requirements of an IMF, h1(t) is denoted as the ith IMF ci(t) and substitutes the residue r1(t) for the original time series data x(t); that is,
(3)r1(t)=x(t)-h1(t).
Otherwise, h1(t) is not an IMF. Then, it substitutes h1(t) for the original time series x(t).Step 6.
Repeat from Step1 to Step 5. The sifting process stops when the residue satisfies one of the termination criteria. First, the residue or the ith component is smaller than the predetermined threshold or becomes a monotonic function such that no more IMF can be extracted. Second, the number of zero crossings and extrema is the same as that of the successive sifting step.By using the above algorithm, the original time series datax(t) can be decomposed into n modes and a residue as follows:
(4)x(t)=∑i=1nci(t)+rn(t),
where n is the number of IMFs, ci(t) represents IMFs which are nearly orthogonal to each other and periodic, and rn(t) is the final residue which is a constant or a trend. By the sifting process, each IMF is independent and specific for expressing the local characteristics of the original time series data. The set of IMFs is derived from high frequency to low frequency, while rn(t) represents the central tendency of data series x(t). In addition, EMD can also be taken as a filter of high pass, band pass, or low pass.


## 2.2. Phase Space Reconstruction (PSR)

Takens embedding theorem [26] provides theoretical foundation for the analysis of time series generated by nonlinear dynamical systems. Later, Sauer et al. [27] show that a phase space can be reconstructed from a univariate chaotic time series. Let a univariate time series {xk}k=1n, where N is the length of the time series, generated from a d-dimension chaotic attractor, then a phase space Rd of the attractor can be reconstructed by using the delay coordinate defined as
(5)x(k)=[xk-τ,…,xk-mτ],
where m is known as the embedding dimension of reconstructed phase space and τ is the delay constant.The selection of the embedding dimensionm and the delay constant τ is very important for prediction modeling [3]. Therefore, for a given delay time τ, a time series {xk}k=1n is represented in the so-called “phase space” by a set of delay vectors (DVs) x(k)=[xk-τ,…,xk-mτ] of a given embedding dimension m.(
1) Determination of the Delay Time
τ
and the Minimum Embedding Dimension
m. The entropy ratio (ER) method [28] is a novel method for determining the set of parameters for a phase space representation of a time series. Based upon the differential entropy, both the optimal embedding dimension and time lag are simultaneously determined.Based upon the probability density function of data, the differential entropy is defined as(6)H(x)=-∫-∞+∞p(x)ln⁡p(x)dx.Particularly convenient is the Kozachenko-Leonenko (K-L) estimate of the differential entropy(7)H(x)=∑j=1Nln⁡(Nρj)+ln⁡2+CE,
where N is the number of samples in the dataset, ρj is the Euclidean distance of the jth delay vector to its nearest neighbour, and CE is the Euler constant. For a given embedding dimension m and time lag τ, let H(x,m,τ) denote the differential entropies estimated for time delay embedded versions of a time series x.The K-L estimates for the time delay embedded versions of the original time seriesH(x,m,τ) and H(xs,i,m,τ) are computed using (7) for increasing m and τ (index i refers to the ith surrogate). To determine the optimal embedding parameters, the ratio
(8)I(m,τ)=H(x,m,τ)〈H(xs,i,m,τ)〉i
needs to be minimized, where 〈·〉 denotes the average over i. To penalize for higher embedding dimensions, the minimum description length (MDL) method is superimposed, yielding the “entropy ratio” (ER)
(9)Rent(m,τ)=I(m,τ)+mln⁡NN,
where N is the number of delay vectors, which is kept constant for all values of m and τ under consideration. The minimum of the plot of the entropy ratio yields the optimal set of embedding parameters.(
2) Identification of Nonlinearity. The delay vector variance (DVV) method [29] is a novel analysis of a time series which examines a signal’s unpredictability by observing the variability of the targets belonging to sets of similar delay vectors (DVs). The DVV method can be summarized as follows for a given embedding dimension.(i)
The meanμd and standard deviation σd are computed over all pairwise distances between DVs.(ii)
The setsΩk are generated, which consist of all DVs that lie closer to x(k) than a certain distance. The distances are taken from the interval [μd-ndσd,μd+ndσd], for example, uniformly spaced, where nd is a parameter controlling the span over which to perform the DVV analysis.(iii)
For every setΩk, the variance of the corresponding targets, σk2, is computed. The average over all sets, divided by the variance of the time series σx2, yields the measure of unpredictability σ*2:
(10)σ*2=(1/N)∑k=1Nσk2σx2.In the following study, the linear or nonlinear nature of the time series is examined by performing DVV analyses on both the original and a number of surrogate time series, using the optimal embedding dimension of the original time series.(
3) Identification of Chaotic Characteristic. When the largest Lyapunov exponent of the system is larger than zero, it indicates that there is a chaotic attractor which can be used to measure the chaotic degree [30]. The largest Lyapunov exponent is computed by the Wolf method [31].


## 2.3. Extreme Learning Machine (ELM)

Suppose there areN distinct samples {(xk,tk)}k=1N, where xk∈Rn and tk∈Rm. The SLFN with M hidden neurons can be described as
(11)∑i=1Mβig(wi·xk+bi)=ok,k=1,2,…,N,
where wi=[ωi1,ωi2,…,ωin]T  is the weight vector connecting the ith hidden neuron and the input neurons, βi=[βi1,βi2,…,βim]T  is the weight vector connecting the ith hidden neuron and the output neurons, oi=[oi1,oi2,…,oim]T  is the actual output vector, and bi is the threshold of the ith hidden neuron. g(·) represents the activation function of hidden neuron, and the operation wi·xk denotes the inner product of wi and xk.If the SLFN can approximate theN samples with a zero error then we have
(12)∑k=1N∥ok-tk∥=0.Thus, there also exist parametersβi, wi, and bi such that
(13)∑i=1Mβig(wi·xk+bi)=tk,k=1,2,…,N.The above equations can be compactly described as(14)Hβ=T,
where
(15)H(w1,w2,…,wM,b1,b2,…,bM,x1,x2,…,xN)=[g(w1·x1+b1)⋯g(wM·x1+bM)⋮⋮⋮g(w1·xN+b1)⋯g(wM·xN+bM)]N×M,β=[β1T⋮βMT]M×m,T=[t1T⋮tNT]N×m.Unlike the traditional function approximation theories which require the adjustment of input weights and hidden layer biases, the input weights and hidden biases are randomly generated. Thus, training an SLFN is simply equivalent to finding a least-squares solution of the linear functionHβ=T. Consider
(16)∥Hβ^-T∥=min⁡β∥Hβ-T∥.
The smallest norm least-squares solution of the above linear system is
(17)β^=H+T,
where H+ is the Moore-Penrose generalized inverse of matrix H. Owing to the Moore-Penrose generalized inverse, the learning speed is dramatically increased for the single hidden-layer feedforward neural network [12].


## 3. Uranium Resource Price Forecasting Method Based on EMD-PSR-ELM

The proposed hybrid approach for international uranium resource price forecasting, namely, EMD-PSR-ELM, combines EMD, PSR, and ELM, and it is composed of four main stages. These four stages are described as follows.Stage 1 (EMD decomposition).
The original time seriesx(t), t=1,2,…,N is decomposed into n IMF components, cj(t), j=1,2,…,n, and one residual component rn(t) by using EMD.Stage 2 (combination of the decomposition components).
In this stage, owing to each IMF with different time scales, high frequency and low frequency components can be obtained by combining IMFs according to the frequency from high to low, and the residue is treated separately. The process can be performed as follows.Step
1. The mean of each IMF is evaluated in order.Step
2. Determining the first IMF ci(t) whose mean significant deviation from zero.Step
3. The first i-1 IMFs are added and reconstructed to the high frequency component c~1(t), namely, c~1(t)=∑k=1i-1ck(t), and the accumulation of the remaining n-i+1 IMFs is reconstructed to the low frequency component c~2(t), namely, c~2(t)=∑k=inck(t).Stage 3 (ELM modeling).
This stage can be subdivided into three steps as follows.Step
1 (phase space reconstruction). Firstly, some parameters such as the delay time and the embedding dimension should be determined. And then, one dimensional time series dataset obtained by combination according to Stage 2 can carry out phase space reconstruction. Thus, the high-dimensional datasets D-={X(t),Y(t)}, t=1,2,…,M can be constructed, where X(t)=[x(t),x(t+τ),…,x(t+(m-1)τ)], Y(t)=x(t+1+(m-1)τ), and t=1,2,…,M, M=N-1-(m-1)τ. Therefore, the input and output samples can be, respectively, represented by the matrix X and Y, in the forms:
(18)X=[x(1)x(1+τ)⋯x(1+(m-1)τ)x(2)x(2+τ)⋯x(2+(m-1)τ)⋮⋮⋮⋮x(M)x(M+τ)⋯x(M+(m-1)τ)],Y=[x(2+(m-1)τ)x(3+(m-1)τ)⋮x(M+1+(m-1)τ)].Step
2 (data normalization). The data needs to be represented in a normalized scale for ELM training and prediction. Thus, in this study, the dataset of three phase space domain are linearly scaled in the range of [0,1] using the following expression:
(19)x-ij=xij-xjmin⁡xjmax⁡-xjmin⁡,
where x-ij represents the jth dimension of ith sample after normalization; xij is jth dimension of ith sample before normalization; xjmin⁡ and xjmax⁡ represent minimum and maximum values of jth dimension, respectively. Then the data was divided into two sets, training set and testing data.Step
3 (ELM prediction). Regression forecast model is set up for the high frequency component, the low frequency component, and the residue by using ELM, respectively.Stage 4 (result composition).
The final prediction results were obtained by compositing the prediction values after denormalization.The proposed EMD-PSR-ELM method is schematically depicted in Figure1.Figure 1
Schematic representation of the proposed EMD-PSR-ELM method.


## 4. Experimental Results

### 4.1. Data Description and Evaluation Criteria

In this paper, for evaluating the performance of proposed EMD-PSR-ELM prediction model, the real time series about international uranium resource prices data were chosen as experimental samples. The data used in this study are monthly data which are freely available from the IndexMundi website (http://www.indexmundi.com) and cover the period from October 1982 to September 2012 with a total of 360 values. Firstly, the time series are analyzed by chaos theory. Delay time τ=4 and embedding dimension m=3 can be simultaneously determined (Figure 2). Secondly, nonlinearity is analyzed by using DVV method. Due to the standardization of the distance axis, these plots can be conveniently combined in a scatter diagram, where the horizontal axis corresponds to the DVV-plot of the original time series and the vertical to that of the surrogate time series. If the surrogate time series yield DVV-plots similar to the original, the “DVV scatter diagram” coincides with the bisector line, and the original time series is probably linear. The deviation from the bisector line is, thus, a measure of nonlinearity; they can be seen from Figures 3 and 4. Additionally, Wolf et al. method [31] is employed to compute the largest Lyapunov exponent λ=0.0339. Therefore, we can conclude that the international uranium resource prices times series has chaotic characteristic owing to λ>0. When these parameters are gotten, the phase space can be reconstructed; that is to say, these optimal embedding dimensions and delay times are used to construct the input matrix. It is easy to know that there are in total 351 data points in the phase space. The data was divided into two sets, training set and testing data, and the first 326 data points are used as the training samples while the remaining 25 data points are used as the testing samples.Figure 2
Plots of the ER for uranium resource price monthly data series.Figure 3
DVV-plots.Figure 4
DVV scatter diagrams.The prediction performance is evaluated using measures including mean absolute error (MAE), root mean squared error (RMSE), and directional prediction statistics (DS) [32]. These measures are as follows:
(20)MAE=1n∑t=1n|x^t-xt|,RMSE=1n∑t=1n(x^t-xt)2,DS=1n∑t=1nat×100%,
where xt and x^t are the actual and prediction values, respectively, n is the sample size, and
(21)at={1,(xt-xt-1)(x^t-xt-1)≥0,0,otherwise.
Obviously, the smaller RMSE, MAE, and larger DS mean better performance.


### 4.2. Forecasting Results

According to previous steps shown in Section3, we carried out the prediction experiments. First, using the EMD technique, the international uranium resource price series can be decomposed into seven independent IMFs and one residue. Figure 5 shows the decomposition results for uranium resource price series using EMD.The IMFs and residue for uranium resource price monthly data series.
(a)
IMF1
(b)
IMF2
(c)
IMF3
(d)
IMF4
(e)
IMF5
(f)
IMF6
(g)
IMF7
(h)
ResidueThrough the analysis of the mean value of each IMF after the decomposition of international uranium resource prices time series data by using EMD, it can be seen from Figure6 that the IMF whose mean significant deviation from zero is IMF3. Therefore, the partial reconstruction with IMF1 and IMF2 represents high frequency components of international uranium resource prices time series data, with the characteristics of small amplitudes, which contains the effects of markets’ short-term fluctuations; and the partial reconstruction with IMF3, IMF4, IMF5, IMF6, and IMF7 represents low frequency components, which should be representative of the effect of these events [33]. The residue is treated as the long-term trend during the evolution of uranium resource prices separately. Figure 7 shows the three components of uranium resource price monthly data series from Oct. 1982 to Sept. 2012.Figure 6
The mean of each IMF of uranium resource price monthly data series.Figure 7
Three components of uranium resource price monthly data series from Oct. 1982 to Sept. 2012.After reconstructing phase space, high frequency, low frequency, and residue are individually used for building ELM prediction models, and the predicted values for the next 25 months uranium resource prices are then obtained by different models. In building EMD-PSR-ELM, the number of hidden nodes is set to 9 for ELM. The error results of the ELM model for the three components are shown in Table1. As can be seen in Table 1, ELM prediction models perform well for all three components.Table 1
The error statistics of the ELM model for the three components.
Component
High frequency
Low frequency
Residue
RMSE
1.3845
0.4799
4.6901
×
10
-
4
MAE
1.0524
0.4181
3.2714
×
10
-
4
DS
64
84
96The forecasting results of the proposed EMD-PSR-ELM model are compared with the PSR-ELM model and PSR-RBFNN model, which uses non-EMD forecasting variables. The PSR-RBFNN forecasting model with three input nodes and one output node is built by using the same phase space reconstruction method as above. The neural network toolbox of MATLAB software is adapted in this study. Mean squared error is 0.1, the spread of radial basis functions is 1.6, and the default settings of neural network toolbox are used for the remaining parameters. In PSR-ELM, the number of hidden nodes is set the same as the above.Figure8 depicts the actual prices and the predicted values from the EMD-PSR-ELM, PSR-RBFNN, and PSR-ELM models. From this figure, it can be observed that there is a smaller deviation between the actual and predicted values using the proposed EMD-PSR-ELM model.Figure 8
Out-of-sample forecasting by different models.Table2 compares the prediction results obtained with the EMD-PSR-ELM, PSR-RBFNN, and PSR-ELM models for the next 25 months uranium resource price.Table 2
Performance of the three forecasting methods for international uranium resource prices.
PSR-RBFNN
PSR-ELM
EMD-PSR-ELM
RMSE
3.4425
2.2389
0.8246
MAE
2.4205
1.4750
0.6919
DS
44
68
80It can be observed from Table2 that the proposed EMD-PSR-ELM model provides a better forecasting result than the PSR-RBFNN and PSR-ELM models in terms of RMSE, MAE, and DS.


## 4.1. Data Description and Evaluation Criteria

In this paper, for evaluating the performance of proposed EMD-PSR-ELM prediction model, the real time series about international uranium resource prices data were chosen as experimental samples. The data used in this study are monthly data which are freely available from the IndexMundi website (http://www.indexmundi.com) and cover the period from October 1982 to September 2012 with a total of 360 values. Firstly, the time series are analyzed by chaos theory. Delay time τ=4 and embedding dimension m=3 can be simultaneously determined (Figure 2). Secondly, nonlinearity is analyzed by using DVV method. Due to the standardization of the distance axis, these plots can be conveniently combined in a scatter diagram, where the horizontal axis corresponds to the DVV-plot of the original time series and the vertical to that of the surrogate time series. If the surrogate time series yield DVV-plots similar to the original, the “DVV scatter diagram” coincides with the bisector line, and the original time series is probably linear. The deviation from the bisector line is, thus, a measure of nonlinearity; they can be seen from Figures 3 and 4. Additionally, Wolf et al. method [31] is employed to compute the largest Lyapunov exponent λ=0.0339. Therefore, we can conclude that the international uranium resource prices times series has chaotic characteristic owing to λ>0. When these parameters are gotten, the phase space can be reconstructed; that is to say, these optimal embedding dimensions and delay times are used to construct the input matrix. It is easy to know that there are in total 351 data points in the phase space. The data was divided into two sets, training set and testing data, and the first 326 data points are used as the training samples while the remaining 25 data points are used as the testing samples.Figure 2
Plots of the ER for uranium resource price monthly data series.Figure 3
DVV-plots.Figure 4
DVV scatter diagrams.The prediction performance is evaluated using measures including mean absolute error (MAE), root mean squared error (RMSE), and directional prediction statistics (DS) [32]. These measures are as follows:
(20)MAE=1n∑t=1n|x^t-xt|,RMSE=1n∑t=1n(x^t-xt)2,DS=1n∑t=1nat×100%,
where xt and x^t are the actual and prediction values, respectively, n is the sample size, and
(21)at={1,(xt-xt-1)(x^t-xt-1)≥0,0,otherwise.
Obviously, the smaller RMSE, MAE, and larger DS mean better performance.


## 4.2. Forecasting Results

According to previous steps shown in Section3, we carried out the prediction experiments. First, using the EMD technique, the international uranium resource price series can be decomposed into seven independent IMFs and one residue. Figure 5 shows the decomposition results for uranium resource price series using EMD.The IMFs and residue for uranium resource price monthly data series.
(a)
IMF1
(b)
IMF2
(c)
IMF3
(d)
IMF4
(e)
IMF5
(f)
IMF6
(g)
IMF7
(h)
ResidueThrough the analysis of the mean value of each IMF after the decomposition of international uranium resource prices time series data by using EMD, it can be seen from Figure6 that the IMF whose mean significant deviation from zero is IMF3. Therefore, the partial reconstruction with IMF1 and IMF2 represents high frequency components of international uranium resource prices time series data, with the characteristics of small amplitudes, which contains the effects of markets’ short-term fluctuations; and the partial reconstruction with IMF3, IMF4, IMF5, IMF6, and IMF7 represents low frequency components, which should be representative of the effect of these events [33]. The residue is treated as the long-term trend during the evolution of uranium resource prices separately. Figure 7 shows the three components of uranium resource price monthly data series from Oct. 1982 to Sept. 2012.Figure 6
The mean of each IMF of uranium resource price monthly data series.Figure 7
Three components of uranium resource price monthly data series from Oct. 1982 to Sept. 2012.After reconstructing phase space, high frequency, low frequency, and residue are individually used for building ELM prediction models, and the predicted values for the next 25 months uranium resource prices are then obtained by different models. In building EMD-PSR-ELM, the number of hidden nodes is set to 9 for ELM. The error results of the ELM model for the three components are shown in Table1. As can be seen in Table 1, ELM prediction models perform well for all three components.Table 1
The error statistics of the ELM model for the three components.
Component
High frequency
Low frequency
Residue
RMSE
1.3845
0.4799
4.6901
×
10
-
4
MAE
1.0524
0.4181
3.2714
×
10
-
4
DS
64
84
96The forecasting results of the proposed EMD-PSR-ELM model are compared with the PSR-ELM model and PSR-RBFNN model, which uses non-EMD forecasting variables. The PSR-RBFNN forecasting model with three input nodes and one output node is built by using the same phase space reconstruction method as above. The neural network toolbox of MATLAB software is adapted in this study. Mean squared error is 0.1, the spread of radial basis functions is 1.6, and the default settings of neural network toolbox are used for the remaining parameters. In PSR-ELM, the number of hidden nodes is set the same as the above.Figure8 depicts the actual prices and the predicted values from the EMD-PSR-ELM, PSR-RBFNN, and PSR-ELM models. From this figure, it can be observed that there is a smaller deviation between the actual and predicted values using the proposed EMD-PSR-ELM model.Figure 8
Out-of-sample forecasting by different models.Table2 compares the prediction results obtained with the EMD-PSR-ELM, PSR-RBFNN, and PSR-ELM models for the next 25 months uranium resource price.Table 2
Performance of the three forecasting methods for international uranium resource prices.
PSR-RBFNN
PSR-ELM
EMD-PSR-ELM
RMSE
3.4425
2.2389
0.8246
MAE
2.4205
1.4750
0.6919
DS
44
68
80It can be observed from Table2 that the proposed EMD-PSR-ELM model provides a better forecasting result than the PSR-RBFNN and PSR-ELM models in terms of RMSE, MAE, and DS.


## 5. Conclusions

This study has presented a forecasting model for uranium resource price by integrating EMD, PSR, and ELM. In terms of the experimental results presented in this study, we can draw the following conclusions.(1)
EMD can fully capture the local fluctuations of data and can be used as a preprocessor to decompose the complicated raw data into a finite set of IMFs and a residue, which have simpler frequency components and high correlations.(2)
On the one hand, although empirical mode composition is an important tool for multiscale modeling, it suffers from mode mixing and mode splitting. To rectify this issue, the Ensemble EMD [34] can be used in the future. In addition, It is also worth trying to employ Multivariate EMD and in particular the Noise Assisted MEMD in order to calculate EMD free of any artifacts [35–37]. On the other hand, other multiscale models, such as the synchrosqueezed transform in the forecasting task, could also be considered [38].(3)
The network topology of the model has an important influence on prediction performance for ELM and RBF neural networks. It is more objective to identify the chaotic characteristic of uranium resource price series and determine the embedding dimension of the reconstructed phase space by quantitative calculation, and then the determined embedding dimension can be served as the numbers of nodes in input layer for the single hidden-layer feedforward neural network and RBF networks.(4)
RMSE, MAE, and DS are used to measure the forecasting accuracy of the forecasting models. The experimental results reveal that the proposed hybrid EMD-PSR-ELM approach outperforms the other models such as PSR-RBFNN and PSR-ELM. Therefore, the proposed method is very suitable for prediction with nonlinear, nonstationary, and highly complex data and is an efficient method for uranium resource price prediction.



---

*Source: 390579-2014-02-20.xml*