Two-Scale Network Dynamic Model for Energy Commodity Processes

In this work, we examine the relationship between different energy commodity spot prices. To do this, multivariate stochastic models with and without external random interventions describing the price of energy commodities are developed. Random intervention process is described by a continuous jump process. The developed mathematical model is utilized to examine the relationship between energy commodity prices. The time-varying parameters in the stochastic model are estimated using the recently developed parameter identification technique called local lagged adapted generalized method of moment (LLGMM). The LLGMM method provides an iterative scheme for updating statistic coefficients in a system of generalized method of moment/observation equations. The usefulness of the LLGMM approach is illustrated by applying to energy commodity data sets for state and parameter estimation problems. Moreover, the forecasting and confidence interval problems are also investigated (U.S. Patent Pending for the LLGMM method described in this manuscript).


Introduction
Understanding the economy evolution in response to structural changes in the energy commodity network system is important to professional economists. The relationship between the different energy sources and their uses provide insights into many important energy issues. The quantitative behavior of energy commodities in which the trend in price of one commodity coincides with the trend in prices of other commodities has always raised the question of whether there is any relationship between prices of energy commodities. If there is any relationship, then what comes to mind is the extent to which one commodity influences the other. Petroleum, natural gas, coal, nuclear fuel, and renewable energy are termed as primary energy components of the energy goods network system because other sources of energy depend on them. Natural gas is usually found near petroleum. Hence, natural gas and crude oil are rivals in production and substitutes in consumption, and energy theory suggests that the two prices should be related. The electric power sector uses primary energy such as coal to generate electricity, which makes electricity a secondary rather than a primary energy source. According to the US Energy Information Administration (EIA), the major energy goods consumed in the United States are petroleum (oil), natural gas, coal, nuclear, and renewable energy. The majority of users are residential and commercial buildings, industry, transportation, and electric power generators. The pattern of fuel usage varies widely by sector [1]. For example, 71% of total petroleum oil provides 93% of the energy used for transportation; 23% of total petroleum oil provides 17% of energy used for residential and commercial use; 5% of total petroleum oil provides 40% of energy used for industrial use; but only 1% of total petroleum oil provides about 1% of the energy used to generate electric power, whereas coal provides 46% of the energy used to generate electric power and natural gas provides 20% of the energy used to generate electric power. This analysis suggests that the strength of interactions between coal and electricity will be stronger than that of crude oil and electricity, or natural gas and electricity.
Energy price forecasts are highly uncertain. We might expect the price of the natural gas and crude oil to follow  Journal of Energy the same trend because they are often found mixed with oil in oil wells and also of the fact that natural gas is often used in petroleum refining and exploration. Recently, Ramberg and Parsons [2] showed that the cointegration relationship between natural gas and crude oil does not appear to be stable through time. They claimed that though there is cointegration between the two energy prices, there are also recurrent exogeneous factors such as seasonality, episodic heat waves, cold waves, and supply interruption from the hurricane affecting the trends in the prices. Brown and Yücel [3] also observed that the price of natural gas pulled away from oil prices in 2000,2002,2003, and late 2005. Oil prices are influenced by several factors, including some that have mainly short-term impacts and other factors, such as expectations about the future world demand for petroleum, other liquids, and production decisions of the Organization of the Petroleum Exporting Countries (OPEC) [1]. Supply and demand in the world oil market are balanced through responses to price movement with considerable complexity in the evolution of underlying supply and demand expectation process. For petroleum and other liquids, the key determinants of long-term supply and prices can be summarized in four broad categories [1]: the economics of non-OPEC supply, OPEC investment and production decisions, the economics of other liquids supply, and world demand for petroleum and other liquids. An understanding of how changes in the price of one energy commodity are expressed in terms of other energy commodity is needed. This would prove to be useful in predicting price behavior over the long run and further facilitates profit-maximizing processes. To check if there is really indeed a relationship between energy commodities, the need to be able to create a model which explains such commodity price relationship over short-and long-time intervals arises. The relationships between energy commodities have been addressed in [2][3][4][5][6][7][8][9][10][11][12][13][14][15]. The error correction model [2-5, 7, 8]     Journal of Energy is the most commonly used model by authors to describe the relationship between energy commodities. Moreover, Hartley et al. [7] have concluded that the U. S. natural gas and crude oil remain linked in their long-term movements. In addition, it is exhibited that there is strong evidence of a stable relationship between these two energy commodities which are affected by short-run seasonal fluctuations and other factors. The rule of thumb [7] has long been used in the energy industry to relate the natural gas prices to crude oil prices. The rule denoted by the 10-to-1 rule states that the price of natural gas is one-tenth of the price of crude oil prices. Similarly, the 6-to-1 rule states that the price of natural gas is one-sixth of the price of crude oil. It has been examined by Brown and Yücel [3] that these two rules do not perform well when used to assess the relationship between US natural gas price and West Texas Intermediate (WTI) crude oil price for the past 20 + years. Moreover, their errorcorrecting model specifies the relationship between the two commodities. Using this model, they show that when certain factors are taken into account, movements in crude oil prices can shape the price of natural gas. Vezzoli [9] in his work applies an ordinary least squares (OLS) regression on the log of natural gas and crude oil prices. Using this model, he was able to conclude that there is a relationship between natural gas and crude oil prices. Bachmeir and Griffin [4] showed that crude oil, coal, and natural gas in the United States have weak cross-cointegration using the error correction model. Ramberg and Parsons [2] show that any simple formula between natural gas and crude oil prices will leave a portion of the natural gas price unexplained. Furthermore, the relationship between natural gas and crude oil using a vector error correction model [2,3] under the  Figure 2: The graph of interaction coefficients κ 1,1 ðm k , kÞ, κ 1,2 ðm k , kÞ, κ 1,3 ðm k , kÞ, κ 2,1 ðm k , kÞ, κ 2,2 ðm k , kÞ, κ 2,3 ðm k , kÞ, κ 3,1 ðm k , kÞ, κ 3,2 ðm k , kÞ, and κ 3,3 ðm k , kÞ.
cointegration between the two energy commodities and other factors such as recurrent exogeneous factors are presented. Villar and Joutz [10] list some economic factors linking natural gas and crude oil prices, while testing for market integration in the United Kingdom during the time when natural gas was deregulated. Asche et al. [6] have integrated the prices of the energy commodities: natural gas, electricity, and crude oil.
Most of this work is centered around the relationship between prices of energy commodities. We are interested in an interdependence of certain energy commodities. Moreover, a nonlinear multivariate interconnected stochastic model of energy commodities and sources with and without external random intervention processes is developed. Also, parameter and state estimation problems of continuous time nonlinear stochastic dynamic process are motivated to initiate an alternative innovative approach. This led to the introduction of the concept of statistic processes, namely, local sample mean and sample variance which further led to the development of an interconnected discrete-time dynamic system of local statistic processes and its mathematical model discussed in Otunuga et al. [16,17]. This paved the way for extending the LLGMM [16] to a multivariate local lagged adapted generalized method of moments case. The parameters in the multivariate model are estimated using the LLGMM method. These estimated parameters help in Table 2: Estimates δ 1,1 ðm k , kÞ, δ 1,2 ðm k , kÞ, δ 1,3 ðm k , kÞ, δ 2,1 ðm k , kÞ,δ 2,2 ðm k , kÞ, δ 2,3 ðm k , kÞ, δ 3,1 ðm k , kÞ, δ 3,2 ðm k , kÞ, and δ 3,3 ðm k , kÞ.  7 Journal of Energy analyzing the short-and long-term relationships between the energy commodities. It has been shown in Appendices A, B, C, and D of the work of Otunuga et al. [16] that the LLGMM parameter estimation scheme performs better than existing nonparametric statistical methods. A performance comparison (with appropriate statistical indices) of the LLGMM method with existing orthogonality condition based on generalized method of moments (OCBGMM) and aggregated generalized method of moments (AGMM) methods using the energy commodity stochastic model is discussed in their work. The method is applied to estimate time-varying parameter estimates in a stochastic differential equation governing energy commodities, stock price processes, and transmission of infectious diseases in the work of Otunuga et al. [16], Otunuga [18], and Mummert and Otunuga [19], respectively. Algorithm for implementing the LLGMM approach is presented in Otunuga et al. [17]. In this work, the usefulness of this approach is further illustrated by applying the technique to study the relationship between three energy commodity data sets: Henry Hub natural gas, crude oil, and coal data sets for state and parameter estimation problems. Interested readers are directed to the work of Otunuga et al. [16][17][18]20] to read more about the LLGMM parameter estimation procedure. Moreover, the forecasting and confidence interval problems are also investigated.
The organization of this paper is as follows.
In Section 2, we derive a multivariate stochastic dynamical model for energy commodities. In Section 3, the multivariate model derived in Section 2 is validated. Due to random intervention, we introduce interventions described by a continuous jump process in our model in Section 4. In Section 5, we discuss the discrete-time dynamic model for sample mean and covariance processes with jump using the work of Otunuga et al. in [16,17]. In Section 6, we discuss about  Figure 3: The graph of interaction coefficients δ 1,1 ðm k , kÞ, δ 1,2 ðm k , kÞ, δ 1,3 ðm k , kÞ, δ 2,1 ðm k , kÞ, δ 2,2 ðm k , kÞ, δ 2,3 ðm k , kÞ, δ 3,1 ðm k , kÞ, δ 3,2 ðm k , kÞ, and δ 3,3 ðm k , kÞ.   10 Journal of Energy parametric estimation techniques. We construct an observation system from a nonlinear stochastic functional differential equation. In Section 7, using the method of moments [21], in the context of lagged adaptive expectation process [22], we briefly outline a procedure to estimate the state parameters for the dynamical model with jump and the model without jump locally. Moreover, the usefulness of computational algorithm is illustrated by applying the procedure to test for the relationship between Henry Hub natural gas, crude oil, and coal for the state and parameter estimation problems. In Section 8, the forecasting and confidence interval problems are also addressed.

Model Derivation
Let p = ½p 1 , p 2 , ⋯, p n Τ be a vector of n energy commodity prices considered to have long-run or short-run relationship with each other, with p j ðtÞ being the price of the jth energy commodity at time t. The economic principles of demand and supply processes suggest that the price of an energy commodity will remain within a given finite expected lower and upper bounds. We define u j ∈ R + = ð0,∞Þ and l j ≥ 0 as the expected upper and lower limits of the jth energy commodity spot prices, respectively. In the absence of interactions between the energy commodities p j , j ∈ Ið1, nÞ, where the market potential for the jth commodity per unit of time at time t can be characterized by ðu j − p j Þðl j + p j Þ. This simple idea leads to the following economic principle regarding the  Figure 4: The graph of interaction coefficients γ 1,1 ðm k , kÞ, γ 1,2 ðm k , kÞ, γ 1,3 ðm k , kÞ, γ 2,1 ðm k , kÞ, γ 2,2 ðm k , kÞ, γ 2,3 ðm k , kÞ, γ 3,1 ðm k , kÞ, γ 3,2 ðm k , kÞ, and γ 3,3 ðm k , kÞ.
dynamic of the price p j of the jth energy commodity. The change in spot price of the jth energy commodity Δp j ðtÞ = p j ðt + ΔtÞ − p j ðtÞ over the interval of length |Δt| is directly proportional to the market potential price, that is, This implies that for some constant α j . From this deterministic mathematical model, if α j > 0, we note that as the price falls below the expected price u j , the price of the jth commodity rises, and as the price lies above u j , there is a tendency for the price to fall. Similar argument follows if α j < 0. Hence u j is the equilibrium state of (3). In a real world situation, the expected upper price limit u j of the jth commodity is not a constant parameter. It varies with time, and moreover it is subjected to random environmental perturbations. Therefore, we consider where ξ j is a white noise process that characterizes the measure of random fluctuation of the upper price limit of the jth commodity; here y j stands for the mean of the energy spot 12 Journal of Energy price process of the jth commodity at time t. It is further assumed that y j is governed by a similar dynamic forces described in (3), that is, where μ j ∈ R + is defined as the mean reversion rate of the mean of the jth commodity. By following the argument used in (4), we incorporate the effects of random environmental perturbations into the lower limit υ j of the mean of the jth commodity: where v j ≥ 0, and e j is a white noise process describing the measure of random influence on the mean price of the jth commodity. Substituting expressions in (4) and (6) into (3) and (5), respectively, we obtain dy j = μ j u j − y j v j + y j dt + μ j u j − y j e j t ð Þdt dp j = α j y j − p j l j + p j dt + α j l j + p j ξ j t ð Þdt: In the absence of interactions and using (7), the system of stochastic model for isolated expected spot and spot prices Table 4: Estimates σ 1,1 ðm k , kÞ, σ 1,2 ðm k , kÞ, σ 1,3 ðm k , kÞ, σ 2,1 ðm k , kÞ, σ 2,2 ðm k , kÞ, σ 2,3 ðm k , kÞ, σ 3,1 ðm k , kÞ, σ 3,2 ðm k , kÞ, and σ 3,3 ðm k , kÞ.  14 Journal of Energy processes are described by the following nonlinear system of stochastic differential equations: where μ j e j t ð Þdt = δ j,j dW j,j t ð Þ, and δ j,j , σ i,j are nonnegative for j ∈ Ið1, nÞ.
Thus, the interconnected energy commodity network system is described by for j ∈ Ið1, nÞ, W j and Z j are n-dimensional independent Wiener processes defined on a filtered probability space (Ω, F, ðF t Þ t≥0 , P ); for l ≠ i, E½dW j,l dW k,i = 0, and for l = i, E½dW j,l dW k,i = dt; the filtration function ðF t Þ t≥0 is right continuous; for each t ≥ 0, each F t contains all P -null sets in F; the n-dimensional random vectors yðt 0 Þ and pðt 0 Þ are F t 0 measurable.
The network system of stochastic differential equations in (13) can be written as follows: and and Yðt, yÞ, σðt, pÞ are a n × n block matrix with each entries having order 1 × n.
In the next section, we outline the model validation problems of (14), namely, the existence and uniqueness of solution process.

Mathematical Model Validation
Here, we validate the mathematical model derived in Section 2. We note that the classical existence and uniqueness theo-rem [23][24][25] is not directly applicable to (14). We need to modify the existence and uniqueness results. The modification is based on Theorem 3.4 [23]. We show the global existence of the solution process of a system of differential equations (14). We note that the rate functions a, b, Y, and σ in stochastic system of differential equations (14) do not satisfy the classical existence and uniqueness conditions [23]. However, these rate functions do satisfy the local Lipschitz condition. Therefore, we construct sequences of functions for the drift and diffusion coefficients of interconnected dynamic system (14) so that the classical conditions for existence and uniqueness theorem are applicable. The construction of modification scheme is as follows: First, we 18 Journal of Energy define a cylindrical subset ½t 0 , ∞Þ × U m of ½0, ∞Þ × R n for t 0 ∈ ½0, ∞Þ, m ≥ 1, where U m is an n-dimensional sphere with radius m defined by The developed stochastic network model (14) can be written as: Here, for each j ∈ Ið1, nÞ and x ∈ R n , we define q j ðx, mÞ Hence, qðx, mÞ = ½q 1 ðx, mÞ, ⋯, q j ðx, mÞ, ⋯, q n ðx, mÞ Τ .
Remark 1. We observe that qðx, mÞ satisfies the global Lipschitz condition on R n with a Lipschitz constant 1.
Using this, together with the local Lipschitz condition assumption on the drift and diffusion coefficients of stochastic differential equations (14), the modified rate coefficient functions in (18) satisfy the classical existence and uniqueness conditions [26,27]. Thus, its solution is denoted by (y m , p m ), for m ≥ 1. Moreover, (y, p) is nonnegative whenever y 0 , p 0 ∈ R n + . Now, we apply Theorems 3.4 and 3.5 of [26] in the context of the modified system of stochastic differential equations (18) and Remark 1 to establish the global existence of solution of stochastic differential equations in (13). For this purpose, we outline the argument used in the proof of these theorems. In addition to the local Lipschitz conditions on the drift and diffusion coefficients, we further impose the following hypothesis on the coefficients: (H 1 ) |ψ j,l t, y ð Þ| ≤ a 2,j +δ j,l y k k, where for i ∈ Ið1, 2Þ, a i,j , a i,j ′ are nonnegative; κ j , γ j ,δ j,l , σ j,l ∈ R + . From (18), we further remark that the dynamic of the mean spot price vector y is decoupled with the dynamic of spot price p. Now, we first apply Theorems 3.4 and 3.5 of [23] in the context of modified system of stochastic differential equations (18) and hypothesis (H 1 ) to establish the global existence of solution of the completely decoupled subsystem of stochastic differential equations in (18). For this purpose, we outline the argument used in the proof of these theorems.

Definition 2.
Let τ m be the first exit time of the solution process y m from the set Bðy 0 , mÞ. Define τ to be the (finite or infinite) limit of the monotone increasing sequence τ m as m ⟶ ∞.
We wish to show that In the following, we present a result that is parallel to Theorem 3.5 [26] in the context of the completely decoupled subsystem of stochastic differential equation (14). For this purpose, it is enough to exhibit the global existence result for the transformed system (18).

Lemma 3.
For m ≥ 1, and y 0 ∈ R n + , let y m ðtÞ = y m ðt, t 0 , y 0 Þ be the solution of the completely decoupled subsystem of (18), and let the hypothesis (H 1 ) be satisfied. Let V 1 be a function defined on ½t 0 ,∞Þ × R n + into R + defined by Then, there exists some constant c 1 > 0 such that where L is the differential operator with respect to (14); e = exp ð1Þ.
Moreover, the global existence of solution of the completely decoupled subsystem of (14) is follows.
Proof. It is obvious that V 1 ∈ C 1,2 on ½t 0 ,∞Þ × R n + → R + . In fact, ð∂V 1 ðt, yÞÞ/∂y j = 2y j /ðy + eÞ, ∂ 2 V 1 ðt, yÞ/∂y 2 j = ð2/ ðy + eÞÞ − ð4y 2 j /ðkyk 2 + eÞ 2 Þ, and ∂ 2 V 1 ðt, yÞ/∂y i y j = −4y i y j / ðkyk 2 + eÞ 2 exist and are continuous functions defined on 19 Journal of Energy ½t 0 ,∞Þ × R n + ⟶ R. Moreover, the L operator with respect to the completely decoupled component is as follows: where Furthermore, inf To show that Probðτ = ∞Þ = 1, we define a function It is obvious that LV ≤ 0. By defining τ m ðtÞ = min ðτ m , tÞ; YðtÞ = y m ðtÞ for t < τ m ; and imitating the argument of Lemma 4 [23], we have Hence, The global existence and uniqueness of the solution of the first component of (18) follows by letting m ⟶ ∞. Hence, from this and the fact that the solution process of transformed system (18) is almost surely identical with the solution process of the original system (14), we conclude the global existence and uniqueness of (14).
Following the idea of Lemma 3, we present a global existence and uniqueness of solution of the system of stochastic differential equations governing the subsystem p in (14).

Lemma 4.
For m ≥ 1, and y 0 , p 0 ∈ R n + , let p m ðtÞ = p m ðt, t 0 , p 0 Þ be the solution of the system of stochastic differential equations governing the subsystem p described in (18), and let the hypothesis (H 1 ) be satisfied. Let V 2 be a nonnegative function on ½t 0 , ∞Þ × R n + into R + defined by 20

Journal of Energy
Then, there exist a constant c > 0 such that where L is the differential operator with respect to (14); e = exp ð1Þ.
Moreover, the global existence of the solution of the system of stochastic differential equations governing the subsystem p described in (14) is as follows.
Proof. It is obvious that and are continuous functions defined on ½t 0 , ∞Þ × R n + ⟶ R. Moreover, the L operator with respect to the system of stochastic differential equations governing the subsystem p described in (14) is obtained as follows: where Furthermore, inf ing the final argument used in proving the global existence of y in Lemma 3, we conclude that there exists a unique global solution to the interconnected system of stochastic differential (14).
In the next section, we discuss the case where we incorporate jump process into the system (y, p).

Energy Commodity Model with and without Jumps
Due to random interventions that affects the price of energy commodities, we introduce random perturbations described by a continuous jump in our model. We follow the approach discussed in [28,29] where a class of stochastic hybrid dynamic processes is investigated.

Journal of Energy
Let K ≥ 0 be the number of jumps on the interval ½t 0 , T, for T > 0. For K ≠ 0, let T 1 , ⋯, T K be the jump times over a time interval ½t 0 , T such that t 0 = T 0 ≤ T 1 <⋯<T K ≤ T, where T i denotes the time at which the ith jump occurred in the system ðy, pÞ. For K = 0, no jump has occurred on the interval ½t 0 , T. We denote the ith subinterval by T i−1 ≤ t < T i . Knowing the global existence and uniqueness solution process of the system (14) on the interval ½t 0 , T, T > 0 in Section 3, for i ∈ Ið1, K * Þ, where K * = K if T K = T, and K * = K + 1 if T K < T and K ≠ 0, we consider the solution process on each subinterval ½T i−1 , T i Þ, between jumps. For i ∈ Ið1, K * Þ and K = 0, we have Ið1, KÞ = ∅ and Ið1, K * Þ = f1g. In this case, we consider the solution process on ½t 0 , T. The interconnected system is governed by the following system of continuous time stochastic process: where ÞÞ is the solution of the system of equation (34); for K ≠ 0 and i ∈ Ið1, K * Þ, Π i and Θ i are jump coefficient matrices. These jump coefficients , nÞ. Following the approach in [28,29], the solution of (34) takes the following representation: Remark 5. In case of no jump, that is K = 0, (34) and (36) reduce to .

Discrete-Time Dynamic Model for Local Sample Mean and Covariance Processes
In this section, we extend the discrete-time dynamic model for local sample mean and variance processes given in Otunuga et al. [16,17] to a multivariate case including the jump process. The development of idea and model of statistic for mean and covariance processes was motivated by the state and parameter estimation problems of the continuous time nonlinear stochastic dynamic model of the energy commodity market network. For this purpose, we need to introduce a few definitions and notations. For more information, we direct the readers to the work of Otunuga et al. [16,17]. For each i ∈ Ið1, K * Þ, let τ i−1 and γ i−1 be finite constant time delays such that 0 < γ i−1 ≤ τ i−1 . If K = 0, then there is no jump and we simply write τ i−1 = τ and γ i−1 = γ. Here, τ i−1 characterizes the influence of the past performance history of state of dynamic process. γ i−1 describes the reaction or response time delays. In general, these time delays are unknown and random variables. These types of delays play an important role in developing mathematical models of continuous time [30] and discrete-time [22,31] dynamic Journal of Energy processes. Based on the practical nature of data collection process, it is essential to either transform these time delays into positive integers or design a suitable data collection schedule or discretization process. For this purpose, for each i ∈ Ið1, K * Þ, we describe the discrete version of time delays of τ i and γ i as follows: Here, ½j:j is the integer part of a real number, and for the sake of simplicity, we assume that 0 < γ i−1 < 1, ð γ i−1 = 1Þ.
For each i ∈ Ið1, K * Þ, the partitions ℙ and ℙ i−1 are defined as follows: In fact, the relationship between a set of jump times fT 1 , T 2 , ⋯, T K g and the partition ℙ defined in (40) is as Using these facts, and noting that if be a finite sequence corresponding to the stochastic process x and partition ℙ i−1 defined in (41). We simply write We also recall the definition of forward time shift operator F [16,17,32]: Moreover, P i−1 k is referred to as the m i−1 k point subpartition of the partition ℙ i−1 in (41) of the closed subinterval Remark 9. We note that for K = 0, that is, no jump, we have where P k is referred as the m k point subpartition of partition ℙ in (40) of the closed subinterval ½t k−m k , t k of ½t 0 − τ, T for k ∈ Ið0, NÞ. (43). This restriction sequence is defined by As m i−1 k varies from 2 to r i−1 + S i−1 + k − 1, the corresponding respective local sequence S m i−1 k ,k at t i−1 k varies from . As a result of this, the sequence defined in (44) is also called a m i−1 k local moving sequence. Furthermore, the average corresponding to the local sequence S m i−1 k ,k in (44) is defined by The average/mean defined in (45)   processes (DTIDMLSMVSP) [16,17] is given by 24 Journal of Energy where d i−1 is the order of the system and The proof of this is given in [16].

Parametric Estimation
In this section, we consider a parameter estimation problem in drift and diffusion coefficients of (34). This is achieved by utilizing the lagged adaptive process [22] and the interconnected discrete-time dynamics of local sample mean and variances statistic processes model (48). For each i ∈ Ið1, K * Þ, we consider a general interconnected hybrid system described by the system of stochastic differential equations: 25 Journal of Energy is the jump coefficient matrix; the jump times T i 's are defined in (34). For each j ∈ Ið1, nÞ, the estimate of the jump coefficient Γ i j is given by Let V ∈ C½½−τ,∞ × R n , R m , and its partial derivatives V t , ∂V/∂x, and ∂ 2 V/∂x 2 exist and are continuous on each interval ½T i−1 , T i . We apply Itô-Doob stochastic differential formula [34] to V, and we obtain where the L operator is defined by For (50) and (51), we present the Euler-type discretization scheme [23]: as the filtration process up to time t i−1 k−1 . With regard to the continuous time dynamic system (50) and its transformed system (51), the more general moments of Δxðt i−1 k Þ are as follows: Journal of Energy , and Τ stands for the transpose of the matrix. From (53) and (54), we have This provides the basis for the development of the concept of lagged adaptive expectation process [22] with respect to continuous time stochastic dynamic system (51). This indeed leads to a formulation of m i−1 k local generalized method of moments at t i−1 k . In the following, we state a result that exhibits the existence of the solution of the system of nonlinear algebraic equations. For the sake of easy reference, we shall state the implicit function theorem without proof.
j,l , and σ i−1 j,l are the system parameters on the jump subinterval ½T i−1 , T i Þ; u i−1 j , κ i−1 j,j , γ i−1 j,j , δ i−1 j,j , and σ i−1 j,j are positive; and for l ≠ j, κ i−1 j,l , γ i−1 j,l ∈ R; δ i−1 j,l , σ i−1 j,l are nonnegative. W and Z are independent standard Wiener process on a filtered probability space (Ω, F, ðF t Þ t≥0 , P) with the properties described in (14). It follows that the interconnected system of stochastic differential equations (56) has 4n 2 + 2n parameters. Also, >0 if y l is cooperating with y j , <0 if y l is competing with y j , = 0 if there is no interaction between y l and y j , j, l ∈ I 1, n ð Þ, ð57Þ where for j, l ∈ Ið1, nÞ, the parameters κ j,l , u j , β j , γ j,l , δ j,l , and σ j,l are the system parameters on the interval ½t 0 , T; u j , κ j,j , γ j,j , δ j,j , and σ j,j are positive; and for l ≠ j, κ j,l , γ j,l ∈ R; δ j,l , σ j,l are nonnegative. For each j ∈ Ið1, nÞ, we pick a Lyapunov function in (51) for (56). Using the Itô differential formula [24], we have Journal of Energy  j ≠ l ∈ I 1, n ð Þ, q ∈ I 1, 2n ð Þ:

Journal of Energy
For each i ∈ Ið1, K * Þ and each j ≠ l ∈ Ið1, nÞ, we define by

Journal of Energy
For every j ∈ Ið1, nÞ, we have We define F 1 = fF 1q g q∈Ið1,n+1Þ , F 2 = fF 2q g q∈Ið1,nÞ , G 1 = fG 1q g q∈Ið1,n+1Þ , and G 2 = fG 2q g q∈Ið1,nÞ . Thus, provided that the determinant of each of the Jacobian matrices J F 1 Þ and J G 2 ðfσ i−1 j,r , σ i−1 l,r g n r=1 Þ are not zero, by the application of Theorem 11 (Implicit Function Theorem), we conclude that for every nonconstant m i−1 k -local sequence , j ∈ Ið1, nÞ, respectively. We illustrate this approach using energy commodities natural gas, crude oil and coal [35][36][37] in the next section.
6.2. Illustration: Application to Energy Commodity. In this subsection, we present an illustration of the derived model on the natural gas, crude oil and coal [35][36][37]. Here, j ∈ Ið1, 3Þ and i ∈ Ið1, K * Þ. Moreover, (56) reduces to Journal of Energy  For each j ∈ Ið1, 3Þ, following the argument used in Section 6.1, we have We also have F 1 = fF 1q g q∈Ið1,4Þ , F 2 = fF 2q g q∈Ið1,3Þ , G 1 = fG 1q g q∈Ið1,4Þ , and G 2 = fG 2q g q∈Ið1,3Þ . For each j ∈ Ið1, 3Þ, it follows that the determinant of the Jacobian of F 1 ðu i−1 j , fκ i−1 j,r g r∈Ið1,3Þ Þ, that is, 3Þ Þ, is not zero provided that all parameters fκ j,r g r∈Ið1,3Þ are not zero or provided the

Computational Algorithm
In this section, we outline computational, data organizational, and simulation schemes. We introduce the ideas of iterative data process and data simulation time schedules in relation to the real-time data observation/collection schedule. For the computational estimation of continuous time sto-chastic dynamic system state and parameters, it is essential to identify an admissible set of local conditional sample average and sample covariance parameters, namely, the size of local conditional sample in the context of a partition of time interval ½T i−1 − τ i−1 , T i . Moreover, the discrete-time dynamic model of conditional sample mean and sample covariance statistic processes in Section 5 and the theoretical parameter estimation scheme in Section 6 motivates to outline a computational scheme in a systematic and coherent manner. A brief conceptual computational scheme and simulation process summary is described below.
7.1. Coordination of Data Observation, Iterative Process, and Simulation Schedules. Without loss of generality, we assume that the real data observation/collection partition schedules   Journal of Energy (41). We present definitions of iterative process and simulation time schedule following the definitions in Otunuga et al. [16].
Definition 14. The iterative process time schedule in relation with the real data collection schedule is defined by where F −r i−1 t i−1 k = t i−1 k−r i−1 is a forward shift operator [32].

Conceptual Computational Parameter Estimation
Scheme. For the conceptual computational dynamic system parameter estimation, we need to introduce a few concepts of local admissible sample/data observation size m i−1 k local admissible conditional finite sequence at t i−1 k ∈ Sℙ i−1 and local finite sequence of parameter estimates at t i−1 k .
43 Journal of Energy for j, l ∈ Ið1, nÞ. This leads to a local finite sequence of parameter estimates at t i−1 k defined on OS i−1 k as follows: The above defined collection is denoted by Thus, y i,s 46 Journal of Energy natural gas, (y 2 , p 2 ) denotes the mean spot and the spot price process of crude oil, and (y 3 , p 3 ) denotes the mean spot and the spot price process of coal. Using the discretized scheme (82), we apply the above conceptual computational algorithm for the real-time data sets namely daily Henry Hub natural gas data set, daily crude oil data set, and daily coal data set for the period of 05/04/2009-01/03/2014, [35][36][37]. Using Δt = 1, ε = 0:001, r = 10, and d = 2, the ε-level suboptimal estimates of the parameters at each real data times are described below for each commodity data sets.
The graph of the diffusion coefficient's parameter for the decoupled dynamical system for y without jump incorporated into the dynamical system are given in Figure 3 below.
Figures 6(a)-6(i) are the graphs of σ 1,1 ðm k , kÞ, σ 1,2 ðm k , kÞ, σ 1,3 ðm k , kÞ, σ 2,1 ðm k , kÞ, σ 2,2 ðm k , kÞ, σ 2,3 ðm k , kÞ, σ 3,1 ðm k , kÞ, σ 3,2 ðm k , kÞ, and σ 1,1 ðm k , kÞ against time t k for the daily Henry Hub natural gas price data set [36], daily crude oil price data set [37], and daily coal price data set, respectively. Table 5 shows the real and simulated estimates for the spot price processes p j ðtÞ, j ∈ Ið1, nÞ. Figure 7 below shows the graph of the real and simulated prices for natural gas, crude oil, and coal data set. Figures 7(a)-7(c) show the graph of the real and simulated spot prices for the daily Henry Hub natural gas data set [36], daily crude oil data set [37], and daily coal data set [35], respectively. The red line represents the real data set pðt k Þ, while the blue line represent the simulated data set p s ðm k , kÞ. Here, we begin by using a starting delay of r = 10. The simulation starts from t r = t 10 . It is clear that the graph fits well. To reduce magnitude of error, we increase the magnitude of time delay. We later compare this result with the case where jump is incorporated into the system. 7.5.2. Relationship between Natural Gas, Crude Oil, and Coal: With Jump Incorporated. In this subsubsection, we analyze the relationship between natural gas, crude oil, and coal with the jump process. Here, we apply the above conceptual computational algorithm described in this section for the realtime data sets, namely, daily Henry Hub natural gas data set, daily crude oil data set, and daily coal data set for the    Figure 13: The graph of σ 1,1 ðm k , kÞ, σ 1,2 ðm k , kÞ, σ 1,3 ðm k , kÞ, σ 2,1 ðm k , kÞ, σ 2,2 ðm k , kÞ, σ 2,3 ðm k , kÞ, σ 3,1 ðm k , kÞ, σ 3,2 ðm k , kÞ, and σ 3,3 ðm k , kÞ for natural gas, crude oil, and coal, respectively.

50
Journal of Energy which an entry of a commodity differ by twice the standard deviation or more from the mean of that commodity. These times are now combined into a single array and sorted out in an increasing order. We give the estimate of the jump coefficient matrices Π i and Θ i defined in (34) in Table 7 below. Table 8 gives the estimates u 1 ðm k , kÞ, κ 1,1 ðm k , kÞ, κ 1,2 ðm k , kÞ, κ 1,3 ðm k , kÞ, u 2 ðm k , kÞ, κ 2,1 ðm k , kÞ, κ 2,2 ðm k , kÞ, κ 2,3 ðm k , kÞ, u 3 ðm k , kÞ, κ 3,1 ðm k , kÞ, κ 3,2 ðm k , kÞ, and κ 3,3 ðm k , kÞ for the decoupled dynamical system for y with jump. Figure 8 shows the parameter estimates u 1 ðm k , kÞ, u 2 ðm k , kÞ, and u 3 ðm k , kÞ for the decoupled dynamical system for y with jump. Figures 8(a)-8(c) are the graphs of u 1 ðm k , kÞ, u 2 ðm k , kÞ, and u 3 ðm k , kÞ against time t k for the daily Henry Hub natural gas price data set [36], daily crude oil price data set [37], and daily coal price data set, respectively. By plotting the real data sets, it is easily seen that the graphs of u 1 ðm k , kÞ, u 2 ðm k , kÞ, and u 3 ðm k , kÞ are similar to the graph of the Henry Hub natural gas, crude oil, and coal data set, respectively. We expect this to happen because u j , j ∈ Ið1, 3Þ, are the equilibrium spot price processes of (5).
as follows:  For the case of no jump, the following graphs show the simulated, forecasted, and 95 percent confidence limit for the daily Henry Hub natural gas data set [36], daily crude oil data set [37], and daily coal data set [35], respectively. Figures 15(a)-15(c) show the graph of the forecast and 95 percent confidence limit for the case where there is no jump for the daily Henry Hub natural gas data set [36], daily crude oil data set [37], and daily coal data set [35], respectively. The figure shows two regions: the simulation region S and the forecast region F. For the simulation region S, we plot the real data set together with the simulated data set as described in Figure 14. For forecast region F, we plot the estimate of the forecast as explained in Section 8.1.1.

Illustration: Prediction/Confidence Interval for Energy
Commodities with Jump. For the jump process, the following graphs show the forecast and 95 percent confidence limit for the daily Henry Hub natural gas data set [36], daily crude oil data set [37], and daily coal data set [35], respectively. Figures 16(a)-16(c) show the graphs of the forecast and 95 percent confidence limit for the case where there is jump for the daily Henry Hub Natural gas data set [36], daily Crude Oil data set [37], and daily Coal data set [35], respectively. Figures 16(a)-16(c) show two regions: the simulation region S and the forecast region F. For the simulation region S, we plot the real data set together with the simulated data set as described in Figure 14. For the forecast region F, we plot the estimate of the forecast as explained in Section 8.1.1.

Conclusion and Future Work
We examine the relationship between the prices of the energy commodities: natural gas, crude oil, and coal. To do this, multivariate stochastic models with and without external random interventions are developed to describe the relationship between the price of energy commodities. The timevarying parameters in the stochastic model are estimated using the LLGMM method. Figures 7 and 14 suggest that the model with jump fits better than the model without jump. For j, l ∈ Ið1, 3Þ, the estimates for the drift interaction parameters γ j,l ðm k , t i−1 k Þ for the case where jump is incorporated suggest that there are definitely interactions between the spot price of these three commodities. As discussed in (57) and (58), the sign of these parameters suggest that there is competition or cooperation between commodities l and j. The estimate of the parameter γ j,l ðm k , t i−1 k Þ in Tables 3 and 10 and Figures 4 and 11 suggests that these commodities either compete or cooperate with each other depending on the time period. We can also describe the relationship between any two commodity j and l, j ≠ l ∈ Ið1, 3Þ, based on the overall average b γ j,l = 1/N∑ N k=1 γ j,l ðm k , t i−1 k Þ. For example, for the case where jump is incorporated, b γ 1,3 = −0:0017 and b γ 3,1 = −0:0095. This suggests that on the average, there is competition between these two commodities. Also, b γ 1,2 = 0:0018 and γ 2,1 = 0:0083. This indicates that on the average, there is cooperation between natural gas and crude oil. Finally, b γ 2,3 = −0:0146 and b γ 3,2 = −0:0013. Therefore, on the average, there is competition between crude oil and coal. In the future, we plan to apply the local lagged adapted generalized method of moments to interconnected nonlinear stochastic dynamic model for log-spot price, expected log-spot price, and volatility process. Also, we plan to incorporate delay in the multivariate interconnected nonlinear stochastic model. We plan to be able to apply the extended local lagged adapted generalized method of moments to other multivariate interconnected nonlinear dynamic model different from energy commodity model.

Data Availability
We will provide the data upon request.

Disclosure
We would also like to mention that this research paper was presented in the Ph.D. Dissertation [15] and at the American Mathematical Society (AMS) conference at San Antonio, Texas.

Conflicts of Interest
The authors declare that they have no conflicts of interest.