Dependence Modeling in Energy Markets using Sibuya-type Copulas

The dependence structure between 756 prices for futures on crude oil and natural gas traded on NYMEX is analyzed using a combination of novel time-series and copula tools. We model the log-returns on each commodity individually by Generalized Autoregressive Score models and account for dependence between them by fitting various copulas to corresponding error terms. Our basic assumption is that the dependence structure may vary over time, but the ratio between the joint distribution of error terms and the product of marginal distributions (e.g., Sibuya’s dependence function) remains the same, being time-invariant. By performing conventional goodness-of-fit tests, we select the best copula, being member of the currently introduced class of Sibuya-type copulas.

The evidence of joint movements in energy markets prices can be possibly explained by the simultaneous impact of increasing economic and political activities in different geographic regions.Traditionally, in a first step of statistical analysis, several single equation time-series models are applied to account for autocorrelation and volatility clustering in marginals.The model variability is huge and depend on many circumstances.
Usually, the return series are modeled by state-space time series (in most of the cases GARCH-type) with suitable marginal distributions for error-terms or by factor models.Following the standard theory, the error terms are assumed to be independent and identically distributed.The empirical analysis shows that this theoretical belief is frequently violated in practice, especially during pre-and post-crisis periods.For example, Alexander (2005) studied dependence between prices for futures on crude oil and natural gas.She concluded that when prices of these commodities are fitted by time-series models, the distribution for the error terms are found to be asymmetric, i.e., the dependence can not be modeled correctly by the bivariate Normal distribution.
In the univariate case non-normality is associated to the skewness and leptokurtosis phenomena, or the fat-tail problem.In the multivariate case, the fat-tail problem can be assigned not only to the marginal distributions but also to the probability of large co-movements of the individual returns.In other words, the tail dependence phenomena is a modern statistical standard to describe the amount of extreme dependence, see Jaschke (2014).
Recently, in a second step, professionals select an appropriate copula family and fit it to filtered log-return series in order to gauge their dynamic interdependence.Since the copula is a function of marginal distributions (applied to error terms), an adequate modeling of individual time series is crucial for estimating the dependence structure between underlying log-returns.For example, the copula approach for modeling energy markets was addressed in Gregoire et al. (2008) who analyzed prices of crude oil and natural gas based on one-month-ahead futures contracts traded on the New York Mercantile Exchange (NYMEX) form July 1, 2003 to July 19, 2006.The authors modeled the log-returns individually as time series, selecting appropriate GARCH-type model for marginal series, and selected the "best" copula to depict the existing dependence between error terms.
As a standard, the researches apply the following copulas in their works: Gaussian copula (without tail dependence), Frank and Student's t-copulas (exhibiting both upper and lower tail dependence), Gumbel and Pareto copulas (having upper tail dependence), Clayton (displaying lower tail dependence), BB1 copula family and few other modifications, see Joe (2015).An implicit assumption in many studies is that the dependence structure between the prices of crude oil and natural gas is constant (static) in time, i.e., the corresponding copula C(u, v) joining marginal log-returns should be time-invariant, meaning that It is well known that the absolutely continuous copula satisfying the last relation is the Clayton's one, see Oakes (2005).It exhibits lower tail dependence (but not upper one).Thus, the Clayton copula can not be used to describe the extreme dependence between markets, for example.Nevertheless (and curiously) there can be found "contributions" whose authors indicate the Clayton copula as the "best choice" for modeling extreme market movements.
Another common misunderstanding can be found in many papers and business reports indicating as a "best" fit Gaussian and/or t-copulas.The Gaussian and t-copulas are popular in applications due to their convenience for computation (being members of the elliptical copula family).The Gaussian copula has no tail dependence, while the t-copula allows different degrees of symmetric upper or lower tail dependence which is determined by the estimated degrees of freedom (d f ) parameter.A smaller d f implies higher tail dependence.As the d f goes to infinity, the Student's t-copula converges to the Gaussian copula.The problem is that the static Gaussian and t-copulas are not time invariant and thus, they are not appropriate to describe the existing dynamic dependence between individual log-returns.
If one googles "Dependence modeling in energy markets with copulas", around 120 related works after 2010 will be listed.Many authors conclude that the overall relationship has increased over time, implying growing tail dependency.Owing to these features, the usual static dependence measures (lower and upper tail indexes) should be considered inadequate after 2008 crisis, indicating a preference to local their versions.We will list and comment several recent contributions: • Reboredo (2011) used copula based GARCH model to study the dependence structure between the crude oil benchmark prices in 12 international crude oil markets.It was found that in times of market stress, crude oil prices in each market tend to exhibit comovements with the same intensity.Moreover, crude oil and natural gas log-returns exhibit non-linear dependence; • Alouia et al. ( 2013) analyzed daily closing prices for the Brent crude oil index and MSCI (Morgan Stanley Capital International) stock market indexes of the six Central East European transition economies (Bulgaria, Czech Republic, Hungary, Poland, Romania and Slovenia) over the period from December 1, 2005 to August 20, 2012, totalizing 1753 observations.The authors investigated the dependence structure between oil and stock market returns through time applying the time-varying copula approach, see, e.g., Cherubini et al. (2012).The conclusion is that bias toward evidence of tail dependence is generated when time-varying parameters in the dependence distribution is not allowed; • Jaschke et al. ( 2012) studied the dependence of extreme events in energy markets.In a deep discussion the authors established that the formal consideration of tail dependence indexes may falsely lead to evidence of asymmetric relation between the returns and suggested to use local version based on tail-copula approach, see Schmidt and Stadtmüller (2006).The empirical investigation focused on the dependence between crude oil and natural gas prices based on one-month-ahead futures contracts traded at NYMEX between July 2, 2007 and July 2, 2010; • Beckmann et al. (2016) analyzed oil prices and exchange rates against the dollar and both experienced long swings over recent decades.The authors focused on the evolution of the relationship between oil prices and dollar exchange rates of 12 oil-exporting and oil-importing countries based on a copula approach.Daily data for two 5 year periods between September 2003 and 2013 were used, taking the crash of Lehman Brothers during 2008 as the dividing point.The growing tail dependence shows that extreme events are more likely to occur simultaneously for both series; • Pan (2014) investigated the tail dependence structure between energy market and stock markets returns in the BRIC (Brazil, Russia, India and China) countries over the period from 12 January 2000 to 28 December 2012.The data set included the Bovespa index (Brazil), the RTSI index (Russia), the BSE 30 Sensitivity index (India) and the Shanghai Composite index (China) returns.Pan (2014) showed that the tail dependence increased rather substantially in the financial crisis of 2008.Moreover, the lower tail dependence for all the paired returns is larger than the upper one.Finally, the tail dependence is found to be the strongest for Russia and the weakest for China; • The vine copula model is another powerful tool to analyze the dependence structure in a multivariate setting, consult Kurowicka and Joe (2011).It allows one to define the structure of relationship between the variables by using expert knowledge, concordance of data, or both.Moreover, it can describe the association between the variables through graphical model, or through what is called pair-copulas.We would mention the applied paper by Kiatmanarochy and Sriboonchittaz (2014)  The chance to depict the existing dynamic dependence between commodities in different geographic regions would increase if one starts with a general model as a base.Our choice is the so-called Generalized Autoregressive Score (GAS) models (also known as Dynamic Conditional Score models), recently introduced by Creal et al. (2013) which incorporate many known time-series models as particular cases.We outline main facts about GAS models in Section 2.1.
As we have seen, the time invariance is not a common copula property and depends on the additional appropriate utility criteria applied by the investigator.For example, all copulas corresponding to distributions possessing bivariate lack of memory property are "time-invariant" in the non-aging world (when the distribution of a random vector coincides with the distribution of its residual lifetime vector), see Kulkarni (2006).
As an alternative and complement to the existing methods, we suggest to use Sibuya-type copulas which incorporate many known copula families.Sibuya-class of copulas generate and represent multivariate distributions preserving the Sibuya's dependence function (being the ratio of joint distribution to the product of marginal distributions, e.g., Sibuya (1960)) with respect to the residual lifetime vector.In fact, the dependence structure may vary over time, but Sibuya's dependence function remains in equilibrium, which is a reasonable utility criteria.The class of Sibuya-type copulas is introduced by Pinto and Kolev (2015) and a short description is presented in Section 2.2.
In Section 3 we apply GAS methodology and Sibuya-type copulas to analyze energy markets log-returns for the data used by Gregoire et al. (2008).It happens that the best fit is achieved with the copula of the extended Gumbels's law (see Example 3) connecting the corresponding error terms.We finish the article with a short discussion.

Methodology
We will introduce briefly two basic tools to be used in our proposal: Generalized Autoregressive Score (GAS) models and Sibuya-class of copulas.

GAS Models
In many settings of empirical interest, time variation in a selection of parameters of a model is important for capturing the dynamic behavior.Creal at al. (2013) introduce a new, general framework for building observation driven approach, where parameters are defined as a function of lagged dependent values, exogenous variables, and past observations.Consult for recent developments http://gasmodel.comLet R t denote the dependent variable of interest (log returns, say), f t the time-varying parameter vector, x t a vector of exogenous variables (covariates), all at time t = 1, ..., n.Define R t = {R 1 , ..., R t }, F t = { f 1 , ..., f t } and X t = {x 1 , ..., x t }.
The available information at time t consists of past observations R t−1 , the time-varying parameters in F t−1 and X t .It is assumed that R t is generated by the observation density p(R t−1 , F t−1 , X t ; θ), where θ is a vector of static parameters.
The updating approach for the time-varying parameter f t is given by the Generalized Autoregressive Score model with orders m and q, to be abbreviated GAS (m, q), as follows where ω(θ) is a vector of unknown parameters.Coefficient matrices A i (θ) and B j (θ) in (1) have appropriate dimensions, while s t is a scaled (link) function, which may be treated as "innovation" depending on set {R t−1 , F t−1 ; θ}.
The score vector s t in ( 1) is specified by where S t−1 is the time dependent scaling matrix (approximated by the inverse of Fisher information matrix usually) and The score terms s t in (1) can be interpreted as Gauss-Newton updating step for every new observation R t that becomes available through time.In addition, the parameter f t is amended in the direction of steepest increase of the log-density at time t.
The main advantage of score s t defined by ( 2) and ( 3) is the possibility of its particular driving mechanism selection.There are several intuitive choices for the scaling matrix S t−1 .For example, • When S t−1 is the identity matrix, i.e., S t−1 = I, the recursion (1) captures models such as the autoregressive conditional multinomial model; • Another natural choice is to set S t−1 equal to the (pseudo)-inverse information matrix based on p(R t−1 , F t−1 , X t ; θ), e.g., We mention two basic properties of GAS (m, q) of bivariate continuous distributions only: • The expectation of the score is zero, i.e., E t−1 [∇ t ] = 0.As a result, s t is a martingale difference sequence; Specific choices for scaling the score vector s t transform the GAS model into well-known observation driven models, see complete list in Creal et al. (2013).Two examples are given below.
Example 1 (GARCH model).Consider the basic model where σ t−1 is a time-varying standard deviation.One can check that the GAS(1,1) model with Example 2 (Linear Gaussian state-space model).Let R t be generated by where f t is the time-varying parameter vector, Λ(.) is the link function, and the scalar α t has a linear dynamic specification The static parameter vector θ incorporates additional fixed and unknown coefficients.In the above state-space model δ is a constant and γ is the autoregressive coefficient.
Alternative frameworks for observation driven models within the exponential family of distributions have been suggested in the literature, including the generalized linear autoregressive models, the generalized autoregressive moving average models and the vector multiplicative error models.In contrast to these proposals, GAS models are able to exploit the complete density structure specified by p(R t−1 , F t−1 , X t ; θ), rather than only means and higher moments.
Thus, GAS models can postulate different dynamics for volatilities from fat-tailed distributions, depending on the underlying observation density p(R t−1 , F t−1 , X t ; θ).Extensions of GAS model to asymmetric, long memory, and other more complex dynamics can be considered as well.Therefore, GAS models allow the formulation of a wide range of new models and gives rise to a number of useful observation driven models that have not been investigated before.

Sibuya-type Copulas
Let S X 1 ,X 2 (x 1 , x 2 ) = P(X 1 > x 1 , X 2 > x 2 ) be the joint survival function of non-negative continuous random vector (X 1 , X 2 ).Define the hazard rate components by The Sibuya-type copulas represent a class L(x; a) of bivariate continuous distributions specified by relation where a 0 , a 1 , a 2 ≥ 0, see Pinto and Kolev (2016).
The class L(x; a) can be characterized as follows.
Theorem 1 (Pinto and Kolev, 2016).If the first partial derivatives of S X 1 ,X 2 (x 1 , x 2 ) exist then r 1 (x 1 , x 2 ) + r 2 (x 1 , x 2 ) = a 0 + a 1 x 1 + a 2 x 2 is fulfilled if and only if the corresponding joint survival function can be represented by where S X i (x i ) are the marginal survival functions, i = 1, 2.
The joint survival function S X 1 ,X 2 (x 1 , x 2 ) given by (4) defines a new aging notion.It is valid only for certain marginal distributions satisfying ∂ 2 ∂x 1 x 2 S X 1 ,X 2 (x 1 , x 2 ) ≥ 0. This inequality implies restrictions for the marginal densities f X i (x i ), which determine the parameter space A of the L(x; a) class as follows If a 0 = f X 1 (0) + f X 2 (0) the bivariate distribution is absolutely continuous and if a 0 < f X 1 (0) + f X 2 (0) the corresponding distribution has a singular component along the line Particular cases of the class L(x; a) are: • All distributions possessing bivariate lack of memory property (in particular, Marshall-Olkin's bivariate exponential distribution, being positive quadrant dependent), e.g., Kulkarni (2006): when a 1 = a 2 = 0; exhibiting negative quadrant dependence: when a 0 = 2 and a 1 = a 2 = θ.The survival copula corresponding to ( 5) is given by • Some members of the Generalized Marshall-Olkin distributions introduced by Li and Pellerei (2011) via stochastic representation where T i 's are independent random variables.Put a It happens that the class L(x; a) can be also characterized by "non-aging Sibuya dependence function" Ω X 1 ,X 2 (x 1 , x 2 ) = S X 1 ,X 2 (x 1 , x 2 ) S X 1 (x 1 )S X 2 (x 2 ) , x 1 , x 2 ≥ 0, being time-invariant with respect to the residual life time vector i.e., Ω X 1 ,X 2 (x 1 , x 2 ) = Ω X 1t ,X 2t (x 1 , x 2 ) for all t > 0.
The class L(x; a) is huge.It contains bivariate distributions that can be symmetric or asymmetric, being absolutely continuous or having a singular component, displaying positive or negative quadrant dependence.Hence, the Sibuyacopula class corresponding to (4) is broad.We will give only an example of a new bivariate distribution extending the Gumbel's distribution (5).
who used the GARCH model and the D-vine copula model to analyze the relationship between the crude oil prices of three different continents: Light crude futures of the NYMEX for North America, Brent crude futures of the Intercontinental Exchange (ICE) for Europe, and Oman crude futures of the Dubai Mercantile Exchange (DME) for Asia.The daily closing prices during the period from 26 December 2008 to 28 June 2013 were used to conduct the analysis.