Forecasting Energy Commodity Prices Using Neural Networks

A new machine learning approach for price modeling is proposed. The use of neural networks as an advanced signal processing tool may be successfully used to model and forecast energy commodity prices, such as crude oil, coal, natural gas, and electricity prices. Energy commodities have shown explosive growth in the last decade. They have become a new asset class used also for investment purposes. This creates a huge demand for better modeling as what occurred in the stock markets in the 1970s. Their price behavior presents unique features causing complex dynamics whose prediction is regarded as a challenging task. The use of a Mixture of Gaussian neural network may provide significant improvements with respect to other well-known models. We propose a computationally efficient learning of this neural network using the maximum likelihood estimation approach to calibrate the parameters. The optimal model is identified using a hierarchical constructive procedure that progressively increases the model complexity. Extensive computer simulations validate the proposed approach and provide an accurate description of commodities prices dynamics.


Introduction
Energy is a principal factor of production in all aspects of every economy.Energy price dynamics are affected by complex risk factors, such as political events, extreme weather conditions, and financial market behavior.Crude oil is a key and highly transportable component for the economic development and growth of industrialized and developing countries, where it is refined into the many petroleum products we consume.Over the last 10 years, the global demand for crude oil and gas has increased largely due to the rapidly increasing demands of non-OECD countries, especially China 1 .Local gas and coal are mainly used in the necessary to setup the analytical models widely used so far.For example, using the historical prices is useful as most of the studies have done ; however, it is by no means enough.In the case of electricity prices, neural networks have been applied to provide a short-term forecast of the System Marginal Price.Szkuta used the three layered Multilayer Perceptron MLP model with backpropagation 10 .The training and testing was performed on deregulated Victorian prices given that in 1997 the market turnover was very large.In addition short-term forecasts drive the analysts decisions and the reduction/increase of generation or demand.
In this paper, we propose an alternative approach to describe the dynamics of energy commodity prices based on machine learning and signal processing, adopting neural network techniques already proven successful in a wide range of forecasting problems 26-31 .More precisely, we apply the Mixture of Gaussian MoG neural networks to estimate, to a given degree of approximation, any probability density.This is possible given that MoG use the Gaussian Mixture Models GMM , Probabilistic Principal Component Analysis PPCA , and latent variables for the singular value decomposition SVD of signal subspaces 32-38 , within a general data-driven estimation of probability density functions for clustering, classification, regression, and the related problems 39-46 .The MoG paradigm adopts a learning process, as for any other neural network model.For instance, in the case of regression, the model parameters are obtained using a training set of input-output samples of the function to be approximated 37, 47-49 .To obtain a good generalization capability of the resulting neural network, a good approximation is mandatory on the samples never involved in the training process.A suitable number of mixture components has to be defined; however, the determination of the optimal number is a very critical problem to be solved, especially when dealing with risk sensitive applications such as medical or financial analysis, since the neural network might be easily overfitted in the case of noisy or ill-conditioned data.For this purpose, a hierarchical constructive procedure for the automatic determination of the mixture components is used; it regularizes the network architecture by progressively increasing the model complexity and controlling the computational cost 50, 51 .
In our knowledge, neural network techniques have been applied only to forecast crude oil and electricity prices.In this paper, we want to study the dynamics of the energy price complex, which has shown large unexpected volatility in the last decade.In order to understand the whole picture the entire complex should be studied.In this context, a powerful tool providing accurate price forecasting is needed.Natural gas, coal, and electricity prices, unlike crude oil, present seasonality features that are usually measured using deterministic techniques.In this paper, we aim to forecast short-term price dynamics in order to be able to adequately measure the existing correlations between the various commodities.To this extent the seasonality component of the gas and coal prices will not affect the results.We apply the MoG approach to forecast crude oil, natural gas, electricity, and coal prices using data for both the European and the US market.The proposed system is trained using daily data collected for the last decade.
The paper is organized as follows.In Section 2, the general framework for time series modeling and forecasting is briefly summarized, while in Section 3, the use of neural networks is proposed as nonlinear regression models suited in this regard.In Section 4, the application of MoG neural networks is illustrated and a hierarchical constructive procedure to train MoG networks is proposed.Extensive computer simulations prove the validity of the proposed approach and provide accurate forecasts of the chosen time series.In Section 5, the numerical results obtained on both reference benchmarks and actual time series are reported and, finally, we draw some conclusions in Section 6.

Basic Concepts for Time Series Modeling and Prediction
Prices related to financial assets are commonly collected as time series uniformly sampled at hourly, daily or monthly rate.Let S t be a generic time series of prices, where t denotes the time index.In order to work with stable data sets for modeling, the financial analysis is usually carried out on the return series y t defined as Given the randomness of prices, a return series is conveniently modeled as a discrete time stochastic process, for which we wish to know the conditional probability density function pdf denoted as f y t |I t−1 y t | I t−1 ; θ , where θ is the vector of pdf parameters and I t−1 is the conditioning set associated with all the available information available prior to time t i.e., past observations and estimated models .Although any pdf depends in general on a set of numerical parameters, to simplify the notation, we will explicit in the following this dependence only when necessary.Almost all the prediction techniques aim to estimate the conditional moments of this pdf, which imply explicit dependence on past observations.In fact, the unconditional moments are relevant to the unconditional distribution f y t y t of the return process and they represent the long-term behavior of the time series, assuming no explicit knowledge of the past.We will assume in the following that all the necessary estimations can be realized in a time smaller than the interval between two consecutive observations, in such a way, we can limit our analysis to a "one-step-ahead" prediction, for which the information I t−1 is available for the analysis prior to time t.Otherwise, the prediction should start earlier, by using the information I t−s , s > 1, to perform a prediction at distance s of the sample at time t.
We consider in this work the reference background of econometrics, where an additive model is used for time series: where μ t is a deterministic component, representing the forecast, and t is a random variable, which takes into account the uncertainty of prediction.In fact, t y t − μ t can be considered as the forecast error, or innovation, and it is in itself a random process.Another hypothesis is that f y where the parameters θ {μ t , σ t can be considered as the volatility of the return series at time t and its accurate estimation, with respect to the implied volatility of the actual process, is mandatory to any financial tool.In general, σ 2 t changes over time and this phenomenon is called heteroscedasticity.So, many models tend to capture this behavior as, in particular, the GARCH model considered in this paper, which is widely used in both theoretical and practical applications.A GARCH P, Q model is based on the following estimation of σ 2 t by using past observations 52 :

2.4
A GARCH P, Q model can be considered as an ARMA process in 2 t and it is generally followed by the "stationary condition": by which some assumptions on the unconditional innovation process t are also possible.Namely, 2.5 is a necessary and sufficient condition for which t is a wide-sense stationary WSS process with E t t 0, Cov k , s k , s E k , s k , s 0 for any k / s, and the unconditional variance is finite and converging to

2.6
This is a common assumption when modeling financial time series, where the forecast errors are zero-mean random disturbances that are serially uncorrelated from one period to the next although not independent, evidently.
A GARCH P, Q model is a generalization of the early ARCH model proposed by Engle in 53 and hence, an ARCH Q model coincides with a GARCH 0, Q model.However, specifying the order P , Q of a GARCH model is not easy and it is still an open problem.Consequently, only low orders are usually adopted in most applications.Nevertheless, several extensions of the original GARCH model have been proposed in the past, by specifying different parameterizations to capture serial dependence on the conditional variance.For instance, some of them are integrated GARCH IGARCH , exponential GARCH EGARCH , threshold GARCH TGARCH or GJR , GARCH-in-mean GARCH-M , and so forth.All of them are able to observe some common characteristics of returns series related to energy commodity prices, in particular volatility clustering, leverage effect, and heavier/fat tails, although they remain weak in capturing wild market fluctuations and unanticipated events.
It is well known that volatility is subject to clustering: large shocks, that is, prediction errors increase the volatility and hence large returns in the next steps are more likely.Volatility clustering is a type of heteroscedasticity accounting for some of the excess kurtosis typically observed in financial data.However, the excess kurtosis can also result from non-normal pdf that happen to have fat tails.In this regard, there are possible other choices for the conditional pdf f y t |I t−1 , the most popular is the Student's t-distribution introduced by Bollerslev in 54 .Certain classes of asymmetric GARCH models can also capture the leverage effect, which results in observed returns negatively correlated with changes in volatility.Namely, volatility rises when returns are lower than expected and vice versa.

Advances in Decision Sciences
The time series model should be completed by a suited hypothesis of the conditional mean μ t as well.A general choice can be based on the linear regression model Bearing in mind 2.2 , y t follows a general ARMA process where usual conditions are given on the eigenvalues associated with the characteristic AR and MA polynomials, in order to ensure stationarity and invertibility, respectively.Generally speaking, ARMA processes make the assumption of white Gaussian noise for t ; however, if an ARMA R, M process is coupled with a WSS GARCH P, Q process satisfying 2.5 , then t meets the condition of weak white noise only.By the way, the condition of strong white noise for t is obtained in the particular case of a GARCH 0, 0 or ARCH 0 process, by which σ 2 t is constant over time and hence t would be an independent and identically distributed i.i.d.Gaussian process.The ARMA model gets in this case a classical homoskedastic process.Energy commodity returns are typically modeled as WSS processes, with constant unconditional mean and constant unconditional variance but nonconstant conditional variance.In the following, we will consider a default association of an ARMA R, M process coupled with a GARCH P, Q WSS process, denoted as ARMA R, M -GARCH P, Q , for which the unconditional variance Var y t y t is a computable function of Var t t obtained in 2.6 .
Some generalizations of the ARMA model are also possible for modeling the con-ditional mean μ t ; for example, the ARIMA model, which is stationarized by adding lags of the differenced series and/or lags of the forecast errors, and ARMA with eXogenous inputs ARMAX , where some "exogenous" or independent variables are added as an explanatory regression data.There are also nonlinear variants to such models as the nonlinear ARMA NARMA and nonlinear ARMAX NARMAX .Finally, we outline that the previous ones are all univariate models.A multivariate generalization is possible by using Vector Autoregression VAR models, which are intended to capture the evolution and the interdependencies between multiple time series.As evidenced in the following, the extension to the multivariate analysis of the proposed GMM approach is straightforward.

Generalized Regression Approach by Neural Networks
The general approach to time series modeling and prediction described in the previous section evidences how both conditional mean μ t and volatility σ 2 t can be estimated through a suited regression problem, which can be compactly defined by the following equation for μ t : 3.1 and by the following equation for σ 2 t : where the orders R, M, P , and Q are fixed in advance they are analogous to ARMA and GARCH models ; ω μ t and ω σ t are the parameter vectors of the regression functions h μ and h σ , respectively, which obviously change over time.
We propose in this paper a new approach to solve the regression problem, and it is characterized by three main features as follows.
i A pure data-driven modeling is adopted, which is suited for neural network learning.
ii The parameters of regression functions are determined simultaneously, since 3.1 and 3.2 are coupled through the values of t determined in the previous time.
iii Nonlinear functions are used for h μ and h σ , in particular by means of neural network models.
By following this approach, we are intended to the direct modeling of the time series dynamics and of its volatility, as actually pursued by GARCH models.However, both linear and nonlinear ARMA-GARCH models are global regression methods that do not involve a parametric basis function expansion of the regression models, similarly to spline functions or the MLP neural network.For this reason, they can be affected by the "curse of dimensionality" problem, since their performances dramatically decrease with the increment of the model order because of the increasing sparseness of data in the input space.
We propose the use of clustering for partitioning the data space, so that clusters of points in significant regions of the data space can be linked directly to the basis functions of a nonlinear regression expansion.All the used neural models allow a form of clustering for nonlinear basis function expansion, this is particularly true for the MoG neural networks based on GMM illustrated in Section 4. Unlike classical neural networks applied to financial time series prediction see, e.g., 55, 56 , MoG neural networks are very suited to clustering for basis expansion and, more in general, for time series particularly noisy or leading to a nonconvex or multivalued regression problems.Classical neural networks, being standard regression models, may fail in giving an accurate description of the observed data and providing their statistical distribution features.MoG networks, instead, succeed in modeling the conditional distribution of the return series, in line to what has been previously proposed in 57, 58 .
From a practical point of view, the problem to be solved is the prediction prior to time t of the sample y t of the return series and the prediction of the related volatility.In the paper, these problems are being also referred to, indifferently, as the estimation of the conditional mean μ t and the conditional variance σ 2 t , respectively.The main data base for the financial analysis is the collection in the past of the price series where price S 0 is the first sample ever collected in the past.By using 2.1 , the data base of returns Usually, the analysis at any time t makes use of a training set D t−1 ⊆ I t−1 , which consists of a limited number N T of previous observations.As explained in next section, D t−1 is determined by means of previous models and predictions, using the samples A prediction process is usually applied for a given number N S of time intervals starting at time t T S , that is, for t T S • • • T S N S − 1 and 1 ≤ N T ≤ T S − 1, where N S also represents the horizon of prediction.Two alternatives are possible in this regard: the prediction models are estimated only once, prior to time t T S ; the prediction models are estimated at every time step, by changing consequently the training set at any time.We will consider in the following the second approach, since the former is a generalization of a Advances in Decision Sciences multi-step-ahead prediction for which suited techniques have been proposed in the literature 59, 60 .
The proposed algorithm for the data-driven regression and forecasting can be summarized in the iteration of the following steps.

Initialization
Let t T S and find the initial conditions for k and σ These values can be inferred by using any adequate technique.We used in this regard an ARMA-GARCH model applied to the samples y k from T S − N T to T S − 1.
Step 1.With the current value of t, determine the training set to be used for the model learning.It consists of two matrices Each row of these matrices is an input-output pattern that can be used for learning.In fact, the first M N columns of D μ t and the first P Q columns of D σ t represent the inputs to h μ and h σ , respectively, for every sample of the training set.The last column of both matrices is the expected value to be estimated in correspondence with every pattern.The last row of matrices holds the most recent observation.
Step 2. Determine, at the current time t, the parameters ω μ t of the regression function h μ by using the training matrix D μ t and an appropriate learning algorithm according to the chosen regression model.Similarly, learn the parameters ω σ t of h σ by using D σ t .For example, if an ARMA-GARCH model is used, the parameters can be estimated by maximum Gaussian likelihood 52 .The specific procedure for MoG neural networks is illustrated in the next section.
Step Once the iteration is stopped, we have N S samples of conditional mean forecast , innovation and volatility pertaining to the time interval where prediction is carried out.The performance of prediction can be evaluated by means of suited benchmarks and error measures applied to the obtained results.A useful collection of such measures will be illustrated in Section 5.

Mixture of Gaussian Neural Network Applied to Time Series Analysis
We introduce in the following the architecture of the MoG neural network to be used for regression in 3.1 and 3.2 .It is based on a GMM that is suited to approximate any pdf as well as any generic function.The model parameters are estimated through a maximum likelihood approach and a constructive procedure is adopted in order to find a suitable number of Gaussian components.

The GMM Model for Regression
Let us consider the estimation of the conditional mean first; we should determine the parameters ω μ t of the regression function h μ : R M → .To simplify the notation, let x ∈ R M be a generic column vector in the input space of h μ and let y ∈ be a generic column vector in the output space of h μ .Although the output is a scalar, it is convenient to keep the vector notation also for y.In fact, the regression approach adopted by MoG neural networks can be immediately applied to a multivariate time series analysis and all the successive considerations are independent of the dimension of the input and output spaces where MoG regression is applied.
The peculiarity of the MoG approach is the estimation of the joint pdf of data f x,y x, y , with no distinction between input and output variables.The joint density is successively conditioned, so that the resulting f y|x y | x can be used for approximating any function 61, 62 .The joint pdf f x,y x, y is based on a GMM of C multivariate Gaussian components in the joint input-output space: where m j and K j are, respectively, mean and covariance matrix of the jth multivariate normal distribution, j 1 • • • C, π j is the prior probability of the component, and The conditional pdf f y|x y | x can be readily determined from 4.1 , that is where Advances in Decision Sciences and g j x is the weighted projection into the input space of the jth component, obtained through the marginal densities: The conditional mean and the conditional variance of f y|x y | x can be calculated easily from 4.3 and 4.4 ; when y is a scalar, f y|x y | x will be a univariate conditional pdf and hence

4.7
The mean obtained using 4.6 corresponds to the least-squares estimation of any output y associated with an input pattern x.The said equation defines distinctly the regression model of a MoG neural network, which has some evident similarities with respect to other well-known neural models such as Constrained Topological Mapping CTM , Radial Basis Function RBF , Adaptive Neurofuzzy Inference system ANFIS , and so on.It is a piecewise linear regression model, which is evidently based on a suitable clustering procedure yielding several regions of the input space, defined by g j x , j 1 • • • C, where the input-output mapping can be locally approximated by the linear functions m j y|x obtained in 4.4 .Moreover, by analyzing 4.5 , we notice that C j 1 g j x 1, for any x.

4.8
As a consequence of this constraint, if the basis functions in the input space are well separated and nonoverlapping then, for any x there exists a Gaussian component q such that g q x ∼ 1 and g j x ∼ 0 for j / q.The architecture of the MoG neural network resulting from the determination of the previous regression model, that is, equations from 4.1 to 4.6 , is reported in Figure 1.
Resuming the original problem of time series forecasting, the MoG model can be used to estimate the conditional mean μ t ; namely, it corresponds to a specific implementation of h μ .Looking in particular at 3.1 , 4.4 , and 4.6 , when x x μ t T , we have 4.9 Nevertheless, we remark that a conditional pdf is obtained in 4.3 , which is just an instance of the generic pdf f y t |I t−1 y t | I t−1 ; θ .As a consequence of this consideration, the estimation of MoG parameters can also be considered as a direct estimation of the parameters θ in a GMM conditional pdf of the time series, which yields a nonlinear model 4.9 for the conditional mean and also a nonlinear heteroskedastic model 4.7 for the volatility.By the way, when C 1, there will be only one Gaussian component in the MoG mixture and hence, f y|x y | x will be a simple normal distribution and the MoG regression will be identical to a linear ARMA model with homoskedastic constant volatility σ 2 t K 1 y|x given by 4.4 .In spite of the last considerations, MoG neural networks can also be used for the direct estimation of the volatility in order to complete the heteroskedastic model of time series introduced in 3.2 .Following the computational scheme illustrated in the previous section, we should determine the parameters ω σ t of the function h σ : P Q → .All the previous equations remain valid; in this case, x ∈ P Q will be a generic column vector in the input space of h σ and y ∈ will be associated with the output space of h σ .For instance, the forecast of volatility can be obtained similarly to 4.9

4.10
The use of clustering in MoG network for kernel-based regression is also suited to capture volatility clustering in energy commodity prices.GARCH models are able to model volatility clustering mainly because the model parameters are estimated repeatedly over time.This is also obtained by our proposed approach.However, the standard GARCH model in 2.4 is only a linear one, which unlikely can capture the clusters present inside the training set.Consequently, this makes very critical the choice of the number of past observations to be used for prediction.Nonlinear GARCH models can alleviate this problem, although using a global nonlinear regression.The training algorithms of MoG neural networks and other similar neural models are intended to find such clusters as a primary goal in their optimization routines.The capability of MoG to find volatility clustering and nonlinear phenomena in the time series analysis will be clearly proved in Section 5.
It is also important to mention that an MoG neural network can generalize the GMM mixture model by using latent variables, which allow a parsimonious representation of the mixture model by projecting data in a suitable subspace; consequently, they are widely used in factor analysis, principal component analysis, data coding, and other similar applications 61 .A full equivalence, under basic conditions, has been proved in 63 between the PPCA clustering performed by the MoG network and the reconstruction of signal subspaces based on the SVD 32 .Therefore, this equivalence allows the adaptive determination of clusters identified on data, in order to maximize the quality of the reconstructed signal i.e., the predicted time series .Furthermore, using known results in the SVD framework, the training algorithm may be performed with the automatic estimation of the latent variable subspace or, equivalently, of the signal subspace.

Training of MoG Neural Network
Without loss of generality, let us consider in the following the regression for conditional mean only.A useful way for determining the parameters of the whole Gaussian mixture f x,y x, y , that is, π j , m j and K j , j 1 • • • C, is based on the maximum likelihood ML approach.It concerns with the maximization of the log-likelihood function: For a fixed number C of mixture components, the ML estimation can be pursued by using different optimization schemes as, for example, the expectation-maximization EM algorithm 47 .It consists in the iteration of two steps, that is, the E-step and the M-step.The essential operations in the case of the GMM 4.1 are summarized in the following.i E-step.With the current estimation of the mixture parameters, the posterior probabilities

4.12
ii M-step.With the current posterior probabilities, the parameters of each Gaussian component, j 1 where K j is updated by using the new values of m j .
The algorithm starts with an initial determination of the mixture parameters.It stops when the absolute difference of the log-likelihood values 4.11 calculated in two successive E-steps is less than a given threshold.

Advances in Decision Sciences 13
When training an MoG network, the main problem is related to the local convergence of the EM algorithm and to the correct determination of the number C of Gaussian components.The former problem mainly depends on a good usually random initialization of the mixture parameters; the latter is a well-known problem, which is directly related to the generalization capability of the neural network.In fact, its performance could be inadequate if the training set is overfitted by an excessive number of Gaussian components.A plain solution to these problems could be based on the use of the EM algorithm with different values of C and with different initializations for each value of C. Once a set of different MoG neural networks has been generated, the choice of the best mixture model can be accomplished by relying on the supervised nature of the prediction problem, that is, by using a cost function measuring the overall generalization capability of the network in terms of its prediction error.
The plain optimization approach still suffers from some serious drawbacks, which mainly depend on the number of different initializations performed for each value of C. The lower is the number of initializations, the lower is the probability to obtain a satisfactory mixture after the EM algorithm.Conversely, the higher is the number of initializations, the higher is the computational cost of the optimization procedure.In order to overcome also these problems, we propose the use of a constructive procedure, where C is increased progressively and only one execution of the EM algorithm is needed for every value of C.This procedure eliminates any random initializations of EM algorithm and, consequently, the necessity to optimize different EM solutions for a fixed C. Thus, the computational cost of the whole training procedure is heavily reduced with respect to a plain optimization approach.The algorithm will stop when C reaches a value representing the maximum complexity allowed to the network.In this way, the training algorithm is structurally constrained, so that overfitting due to the presence of outliers in the training set can be prevented when C is low, ensuring robustness with respect to noise in the time series to be predicted.
The constructive procedure is intended to find a component in the GMM performing a poor local approximation and to substitute it by a pair of new components.The underlying idea is to prevent the initialization of the new components in underpopulated zones of the training set, since this is the typical situation where the EM algorithm will converge to a poor local maximum of 4.11 .Several heuristics are possible to select the component to be split; we use in the following the one clearly illustrated in 39 that is based on the supervised nature of the present regression problem.More details about the constructive procedure and the demonstration that it yields better mixture weightings than a random initialization approach can also be found, for instance, in 30 .

Performance Evaluation
The validity of the proposed approach has been validated by extensive computer simulations.Some illustrative examples are summarized in this section, firstly simulated data as reference benchmarks are considered, then actual return series related to energy commodity prices are used.The numerical results are obtained by using well-known neural and neurofuzzy models, which are compared with respect to the commonly used combination of ARMA-GARCH models estimated by maximum Gaussian likelihood.
The training procedure of neural regression models follows the scheme illustrated in Section 3. Bearing in mind the introduced notation, let y t , t T S • • • T S N S − 1 , be the set of actual samples of the time series to be predicted.For each sample, the neural networks are trained by using the previous N T observations of the time series, that is, from t − N T to t − 1, together with the related innovations and conditional variances.A one-step-ahead prediction is therefore applied and this procedure is iterated N S times.After every prediction, the sequences of innovations, conditional mean and conditional variance forecast are updated.
In addition to the suggested MoG neural network, trained by the constructive procedure previously illustrated, we further consider two well-known neural architectures: RBF and ANFIS.The former is a feed-forward neural network trained by a constructive procedure that iteratively creates a radial basis network one neuron at a time.Neurons are added to the network until an error goal or a maximum number of neurons is reached 64, 65 .The ANFIS is a neurofuzzy network, combining the data-driven learning of neural networks with the robustness of fuzzy logic and linguistic fuzzy rules.ANFIS networks have been trained by a subtractive clustering method for rule extraction 66 , while the rule parameters are obtained by means of a standard least-squares method coupled with the back-propagation optimization 67 .All these training procedures also aim to optimize the structural complexity i.e., number of kernels, hidden nodes, fuzzy rules, etc. of the resulting neural network.For reasons of conciseness, we will not provide details about complexity in the following, since it is optimal as long as the neural network exhibits a good generalization capability, which is evaluated as described in the following by means of the network's performance on test sets not used during training.
We take particular care to the criteria used to evaluate the performance of the algorithms.There are many error measures adopted in the literature, in particular for measuring the prediction error.They are differently used according to different fields as statistics, engineering, econometrics, and so forth.Let μ t , t T S where y is the average of the actual N S samples of y t .In general, a prediction error can be defined on the return series only by using the estimated conditional means.However, when dealing with simulated data, we can also determine the true sequence of conditional variances.So, the same measures can be applied to the sequence of predicted conditional variances, considering σ 2 t in place of μ t and the actual conditional variance in place of y t .For this reason, we will distinguish the errors where necessary.
Best approximation is linked to small values of 5.1 , 5.2 , and 5.3 and it is interesting to notice that the MSE is the energy of the prediction error, while in NMSE and NSR, this energy is normalized to the variance and to the energy of the observed sequence, respectively.Each measure aims to quantify how much the prediction has been corrupted by errors and should allow an objective evaluation of the prediction accuracy, independently of the used time series.So, we should be able to set a sound threshold of the NSR that assures a good performance especially for what concerns the statistical distributions.In this circumstance, we will prove that good performances are usually obtained when the NSR is lower than approximately 1 dB.
Energy prices are critical to producers, consumers, and investors.Derivative instruments are an efficient tool for hedging and managing the price and volatility risks.To properly price energy derivatives is critically dependent on accurately modeling the underlying spot price.A successful model describing spot prices dynamics must capture the statistical features of the analyzed data in the simulated series.To this aim, the unconditional moments from the first up to the fourth order are estimated as time averages and considered for both actual and predicted sequences.A given model is suited to forecast and model energy commodity prices when the first four moments of the predicted sequences result as close as possible to the moments estimated on the market data.Being able to reproduce the probability distribution of the observed series together with an accurate prediction of the daily prices will allow investors and risk managers to estimate the profit/loss scenarios to set up the adequate risk management strategies.

Synthetic Benchmarks
For the sake of comparison, we consider the following artificial data sets with Gaussian innovations investigated in 55, 68 , also useful to assess the proposed approach in case of high degree of nonlinearity:

5.7
A realization of 1200 samples is generated by using each model; N T 1000 samples are used to train the prediction model for each sample; N S 200 samples are used for test, starting at sample T S 1001; hence, the index of predicted samples ranges from t 1001 to t 1200.We generated 50 different realizations for both models and we report in the following the mean results over the different realizations.An illustrative realization of mean and variance obtained by the first model is shown in Figure 2, while a realization of the second model is reported in Figure 3.As default model parameters we adopted R 1, M 0, P 1, and Q 1; therefore, the standard reference model is an ARMA 1,0 -GARCH 1,1 , which will be briefly denoted in this context as "GARCH." The numerical results for the two models are summarized in Tables 1 and 2, respectively.As explained, we can apply the error measures to the predicted time series and to the predicted variance or volatility as well.In all cases, when the original time series are highly nonlinear, we find that GARCH models are not able to capture the real nature of the data.In the first benchmark, which is a classic heteroskedastic process, the nonlinear models of neural networks, especially RBF and ANFIS, are able to well represent the sinusoidal mean 5.4 .Although GARCH and MoG obtain comparable results on the mean, MoG is much more accurate to model the volatility since GARCH tends to shape the variability due to mean as an increased volatility 55 .Poor results are obtained also when GARCH is used to estimate the variance of the time series in the second benchmark.The nonlinear regression models provided by neural networks perform better and, in particular, MoG network achieves the best prediction of the variance as illustrated in Figure 4.

Energy Commodity Prices
We consider the time series obtained from the daily prices of specific energy commodities.We carried out several tests to assess the validity of the proposed approach.Coal, natural gas, crude oil, and electricity prices over the period 2001-2010 were collected for both the European and the US markets.The results reported in this paper refer to the US market only, similar results have been obtained for the European data set not showing significant differences in the forecasting ability.The studied commodities and the related indexes are as follows: coal DJUSCL, in $/ton , Henry Hub natural gas HH, in $/MMBtu , crude oil WTI, in $/barrel , and electricity PJM, in $/MWh .For electricity prices, we chose the peak load contract referred to h. 12:00.For each commodity price log-returns are estimated using 2.1 .We chose a well representative time window across the "very critical" year 2008, that is, from the beginning of 2006 to the end of 2009.So, taking into account that we have approximately 250 prices and 249 returns per year in the case of coal, natural gas, and crude oil series, each return series consists of about 1000 samples.In the case of electricity prices we have a series of 365 data, given that electricity is traded every day of the year; in this application, for comparison purposes, we adjust the series to 250 trading days only.Each model is trained on the previous N T 500 samples almost two years and N S 500 samples are predicted, that is, the last two years starting from t 501 up to the last sample of 2009.The prediction errors of the conditional mean are evaluated; in addition, the four unconditional moments are estimated for both the predicted sequences and the related original series.Prior investigations can be made in order to find the best combination of the model orders and the size of the training set as well.The goal of this paper is to propose a new machine learning approach to forecast the time series related to energy commodities, so we only demonstrate how it favorably compares with respect to other choices.Fine tuning for the optimal estimate of any models can be addressed in future research works.However,  in order to obtain an accurate comparison of the performances obtained by the proposed neural networks with respect to standard models, we carried out a preliminary optimization of the main model parameters, that is, R, M, P , and Q, so as to obtain the best performance of the reference "GARCH" model for a given time series.Then, every model will use the same parameters when applied to the same time series.The coal DJUSCL daily prices and returns are shown in Figure 5 and the numerical results are reported in Table 3.The optimal parameters are R 1, M 1, P 5, Q 1; hence, a GARCH reference model ARMA 1, 5 -GARCH 1, 1 is fitted.All the neural models score a prediction error better than GARCH: RBF obtains the best NSR but the skewness is not properly matched as in the case of ANFIS.MoG neural network achieves a good NSR performance of 0.18 dB, and its results are plotted in Figure 6.Although not the best model, MoG neural network is able to follow the dynamic of changes of the actual conditional mean and to forecast the increasing volatility as proved by the behavior of the estimated conditional variance.
The HH daily prices and returns of natural gas are shown in Figure 7 and the numerical results are reported in Table 4.The optimal parameters are in this case R 2, M 2, P 2, Q 1, so a GARCH reference model ARMA 2, 2 -GARCH 2, 1 is fitted.The MoG neural network has the best NSR performance of 0.21 dB and the related moments adequately fit with those of the original time series.A sufficient accuracy is also obtained by RBF and ANFIS neural networks.GARCH is not suitable for the prediction of HH returns, since the moments are estimated very poorly, especially the kurtosis.The numerical results  of the MoG neural network are qualitatively confirmed by the accurate prediction of the conditional mean and volatility reported in the plots of Figure 8.The behaviors of crude oil WTI daily prices and returns are shown in Figure 9.The large volatility at the end of 2008 is the feature that requires an accurate forecasting technique.A more complex model is therefore necessary, using R 4, M 2, P 2, Q 3. The GARCH reference model ARMA 4, 2 -GARCH 2, 3 is evidently outperformed by the neural networks, as evidenced by the results summarized in Table 5.The best NSR is now obtained by ANFIS although the predicted sequence does not fit the skewness of the original one.Once again MoG is able to fit the original moments, also maintaining a good prediction accuracy and following the changes of volatility in the underlying process as illustrated in Figure 10.
Finally, the prices and returns of PJM electricity index are shown in Figure 11 and the numerical results are reported in Table 6.The model parameters are R 3, M 2, P 1, Q 1; the GARCH reference model is ARMA 3, 2 -GARCH 1, 1 .Similarly to the previous cases, MoG neural network performs better than GARCH, RBF, and ANFIS.Globally neural networks improve the NSR performance of more than 2 dB with respect to GARCH despite a biasing that shifts the estimate of the mean to negative values.
PJM returns are not characterized by the same amount of volatility than WTI or coal.However, some power markets usually show some extreme behaviors, in particular spikes that enhance the kurtosis to very high values.Actually, 2008 is deeply affected by these events and it has been analyzed, though compounded with 2009.Anyway, in order to test the proposed approach also in such situations, we focus the last quarter Q4 of 2008 where the kurtosis increases to about 9. Leaving the same the training process, as previously illustrated,       the numerical results obtained in the considered subperiod are summarized in Table 7.These results prove that the proposed neural approach is able to model complex time series showing high kurtosis, as in the presence of spikes, especially when compared to GARCH models.Nonetheless, further comparisons of neural models with respect to stochastic models such as jump-diffusion, regime-switching, etc. can be considered in the future, taking into account the potentiality of neural networks to obtain good performances independently of the considered process, whilst stochastic models must be carefully chosen a priori also relying on the skill of expert users.

Conclusion
In this paper, a new neural network approach is proposed for modeling time series associated with energy commodity prices.It is based on a pure data-driven estimation of nonlinear regression models, which are applied to the estimation of both conditional mean and volatility.In this regard, we focused on MoG neural networks, which are based on GMM for kernel-based regression and pdf approximation.MoG neural networks can be considered as a nonlinear generalization of GARCH models, for which volatility clustering can be pursued in a straightforward manner.The use of a constructive procedure to train MoG neural networks has been proposed.It aims to the regularization of the network complexity, determining automatically the optimal number of components in the GMM in order to avoid overfitting and maximize the generalization capability of the network.The proposed approach provides very accurate description of energy prices dynamics, allowing to estimate daily prices for energy commodities over a long time horizon.The validation performed on historical data shows that the neural network approach generates prices that are able to replicate the daily data and to reproduce the same probability distribution for the various series.
The proposed approach not rely on information about seasonality in the processed time series.Actually, both natural gas and coal exhibit seasonality in the price and volatility, thereby suggesting that there is a nonrandom component to the signals which can be used to improve backtesting of the model using daily returns.However, seasonality is not explicitly accounted for using neural networks, since they can be used as a black-box where such an information is inherently processed using neither specific preprocessing tasks nor the personal skill of the users.Nevertheless, unlike coal, natural gas, and crude oil, electricity is nonstorable and exhibits large jumps in daily spot price returns.Also in this case, the proposed approach and the generalization capability of the considered neural models are able to take into account these phenomena without using specific models as, for example, the well-known mean-reverting jump diffusion model.
The accurate prediction of daily prices on a short-term basis will provide a useful tool to estimate the relationship existing between related commodities, that is, correlation between natural gas and crude oil or between natural gas and electricity.These relationships represent a key element for the pricing of energy derivatives as well as spread options which have become a largely used instrument to efficiently manage risk in the energy markets.Currently, we are investigating more advanced techniques for the application of the proposed approach to a multivariate time series analysis and for the automatic and more reliable selection of the samples to be used for prediction, including the order of regression models and the resulting complexity of neural models.

g 1 Figure 1 :
Figure 1: Regression model obtained by an MoG neural network.The number of scalar inputs, that is, the dimension of the input vector x, is limited to M for illustration.

10 Figure 2 : 5 Figure 3 :
Figure 2: Realization of mean and variance for the sinusoidal mean benchmark defined by 5.4 and 5.5 .

Figure 4 :
Figure 4: Prediction of conditional variance of the zero-mean benchmark obtained by the MoG neural network: forecast black ; actual gray .

2 DJUSCLFigure 5 :
Figure 5: Daily prices and returns of the DJUSCL coal index.

Figure 9 :
Figure 9: Daily prices and returns of WTI crude oil index.

Figure 10 : 2 Figure 11 :
Figure 10: Prediction of crude oil returns using the MoG neural network.
2 t } are determined on the basis of the conditioning information I t−1 , being μ t and σ 2 t mean and variance of the univariate normal distribution.Thus, the conditional mean of y t is E y t |I t−1 y t μ t and the innovation process has a zero-mean normal conditional distribution for which E t |I t−1 t 0 and Var t |I t−1 t E t |I t−1 3. By means of the parameters ω t y t − μ t , t ← t 1, and go back to Step 1 if t ≤ T S .
• • • T S N S −1 , be the set of conditional means representing the prediction obtained by using any model of the corresponding values y t .The error measures used in this paper are the following ones:

Table 1 :
Numerical results for the sinusoidal mean benchmark of 5.4 and 5.5 .

Table 2 :
Numerical results for the zero-mean benchmark of 5.6 and 5.7 .

Table 3 :
Numerical results for DJUSCL return series.

Table 4 :
Numerical results for HH return series.

Table 5 :
Numerical results for WTI return series.

Table 6 :
Numerical results for PJM return series.

Table 7 :
Numerical results for PJM return series in Q4 2008.