A conditionally heteroskedastic model with time-varying coefficients for daily gas spot prices

A novel GARCH(1,1) model, with coefficients function of the realizations of an exogenous process, is considered for the volatility of daily gas prices. A distinctive feature of the model is that it produces non-stationary solutions. The probability properties, and the convergence and asymptotic normality of the Quasi-Maximum Likelihood Estimator (QMLE) have been derived by Regnard and Zakoian (2009). The prediction properties of the model are considered. We derive a strongly consistent estimator of the asymptotic variance of the QMLE. An application to daily gas spot prices from the Zeebruge market is presented. Apart from conditional heteroskedasticity, an empirical finding is the existence of distinct volatility regimes depending on the temperature level.


Introduction
Following the deregulation of natural gas markets in Europe, natural gas supplying historically ran by long term contracts indexed on crude oil between producing countries and retailers was diversified trough new financial markets (National Balancing Point in UK, Zeebrugge market in Belgium), where it can be freely sold at different time horizons. This restructuring has generated uncertainty, requiring the development of appropriate valuation and risk management strategies.
Such strategies require an appropriate modeling of the prices volatility. The standard GARCH models of Engle (1982) and Bollerslev (1986), which arguably constitute the most important class of models for financial data may be inadequate for energy prices. The reason is that energy prices are subject to pronounced daily seasonal patterns, which may not only concern the conditional mean but also the volatility. The periodic ARCH model introduced by Bollerslev and Ghysels (1996) is able to capture those seasonal behaviors in the conditional variance. However, in this model the different regimes appear in a purely periodic succession and it may be worth introducing more flexibility. A GARCH model with regression effects and scaled by seasonal factors has been recently proposed for electricity prices by Koopman, Ooms and Cornaro (2009).
The purpose of this article is to develop a new class of volatility models, introduced in a companion paper by Regnard and Zakoian (2009) (hereafter RZ), for characterizing the seasonal patterns induced by other variables such as temperature. In this model, the parameters associated with the volatility dynamics depend on an exogenous variable, similarly to papers dealing with the conditional mean by Azrak and Mélard (2006), Bibi and Francq (2003), Gautier (2004a, 2004b).
The article is organized as follows. Section 2 introduces the model and its main probability properties. It is shown how the model can be used for prediction purposes and the QML (Quasi-Maximum Likelihood) estimation is discussed. Section 3 proposes an application to gas prices. A preliminary treatment based on a vector error correction model, involving daily gas prices, brent prices and the temperature, is discussed. Finally, the proposed model is fitted with up to five volatility regimes depending on the temperature. The different specifications are tested, and compared via out-of-sample predictions. Section 4 concludes. A technical proof is given in the appendix.

A nonstationary GARCH(1,1) model
The model we consider in this paper is given by where (η t ) is a sequence of independent and identically distributed (iid), centered variables with unit variance; (s t ) is the realization of a process (S t ) with values in a finite set E = {e 1 , . . . , e d }; the functions ω(·), α(·), β(·) are defined on E with values in R + with ω(·) > 0.
In our application, s t will correspond to a level of temperature, observed at time t. For each level of temperature, the volatility is that of a standard GARCH(1,1) model. Thus, if this level remains constant in some period, the volatility is governed by a standard GARCH. When another level of temperature is reached, the specification of the volatility changes. The existence of different regimes for the volatility, is a common feature between this model and the so-called Markov-switching GARCH models (see ...). However, the models interpretations are completely different. In Markov-switching models, the mechanism of regime change is governed by an non observable variable. In our model, it is governed by an observable process which is exogenous to the model. The dynamics of t is conditional to (S t ).

Probability properties
The probabilistic properties of this model have been established by Regnard and Zakoïan (2009). Assuming A0: (s t ) is a realization of a process (S t ) which is stationary, ergodic, defined on the same probability space (Ω, A, P) as (η t ), and independent of (η t ), and letting π j = P (S t = e j ), j = 1, . . . , d and a(x, y) = α(x)y 2 + β(x),

RZ established that if
Model (1) admits a nonanticipative nonexplosive solution ( t ). When γ 0 > 0, the process is explosive: for any initial value σ 2 0 , we have In the ARCH(1) case (no coefficients β), Condition (2) takes the more explicit form d j=1 α π j (j) < e −E log η 2 0 . It can also be noted that the stability of the GARCH(1,1) in each regime, that is is sufficient (but not necessary) for the global stability. A necessary condition for (2) is given by d j=1 β π j (j) < 1. Moreover

Predictions of the squares
For standard GARCH(1,1) models, the optimal prediction of 2 t in the L 2 sense, E( 2 t | { 2 t− , > 0}), is obtained from the ARMA(1,1) representation for the squares. Similarly, for Model (1), This representation is valid because (2) implies from which the existence of the infinite sum in (4) is deduced, by the arguments used to establish the stationarity condition. Note that the expectation of u t conditional on 2 t past values is zero. The optimal predictorˆ 2 t of 2 t , in the L 2 sense, is then Predictions at higher horizons can be derived similarly. Contrary to standard GARCH models, predictions formulas are time dependent trough the coefficients s t .

QML Estimation
The consistency and asymptotic normality of the QMLE have been proven under mild conditions by Berkes, Horváth and Kokoszka (2003) and Francq and Zakoïan (2004) for standard GARCH and ARMA-GARCH models. NZ showed that these properties can be extended to the model of this paper under assumptions which we now detail.
Let θ denote the vector of parameters, with true value θ 0 . The parameter is assumed to belong to a parameter space The sequence (s t ) is observed, and the orders p, q and d are known a priori. Let ( 1 , . . . , n ) be a realization of length n of the nonanticipative solution ( t ). Conditionally on initial values˜ 0 andσ 2 0 the gaussian quasi-likelihood is given by where theσ 2 t are defined recursively, for t ≥ 2, bỹ wherel Indexing the true parameter values by 0, we make the following assumptions.
Then, RZ showed that is a positive-definite matrix, and (σ 2 S,t (θ 0 )) is the process obtained by replacing s t by S t in the second equation of (1).
The following examples illustrate the variability of the asymptotic covariance matrix, for given ARCH(1) coefficients with two regimes, depending on the distributions of (S t ) and (η t ). Let the model and suppose that (S t ) is a Markov chain with transition probabilities p(i, j).
• p(1, 1) = p(2, 2) = 0.95, η t is distributed as a mixing og Gaussian distributions (with κ η ≈ 9): where these asymptotic covariance matrices have been obtained numerically, by simulation of Model (8). The presence of asymptotic covariances equal to zero for parameters of different regimes is due to the absence of coefficients β in the model.

Consistent estimation of the asymptotic variance of the QMLE
To build tests and confidence intervals for the parameters of Model (1), it is essential to have a consistent estimator of the asymptotic covariance matrix of the QMLE. In view of (7), this matrix depends on the distribution of (S t ) which is unknown. However, the following result provides a consistent estimator which can be easily computed.
Proposition 1 Under Assumptions A0-A7, a strongly consistent estimator of the matrix J is given bŷ and a strongly consistent estimator of (κ η − 1)J −1 is Proof. See appendix.

Application to gas prices volatility
We now turn to an example with real data, namely the daily series of gas spot prices from the Zeebruge market. Before modeling the volatility we filter the series from the conditional mean. To capture the joint behavior of the series of gas, Brent prices and temperatures, we consider a VAR model.
We have a sample of daily prices and temperatures from April, 5, 2000 through March, 3, 2008. Figure 1 displays a plot of the gas and Brent series, (G t ) and (B t ) on the whole period. Let g t = log G t and b t = log B t denote the log prices and let (T t ) denote the temperature series. Augmented Dickey-Fuller and KPSS (Kwiatowski, Phillips, Schmidt and Shin, 1992) unit-root tests not reported here suggest that the series g t , b t and T t are integrated of order one.
To filter the gas price conditional mean from the influence of the Brent oil price and the temperature, we use a vector error correction model (VECM). There is a growing literature examining the cointegration relationships between different energy prices. Asche, Osmundsen, and Sandsmark (2006) discuss the cointegration between UK natural gas, Brent oil and electricity prices before and after the opening of the Interconnector in 1998. Bachmeier and Griffin (2006) found evidence of cointegration between crude oil, natural gas and coal in the USA. Panagiotidis and Rutledge (2007) found evidence of a cointegration relationship between the UK wholesale gas prices and the Brent over the period 1996-2003, contradicting the assumption that gas prices and oil prices are decoupled since the liberalisation of gas markets in Europe.

A VECM for gas and brent prices
We begin the analysis with an error correction approach. Recall that, in Johansen's (1988, 1995) notation, a p-dimensional VECM takes the form where ∆ is the difference operator, y t is a p × 1 vector of I(1) variables, µ is a drift parameter, (u t ) is a white noise, Π = αβ is a p × p matrix where α and β are p×r full-rank matrices, with β containing the r cointegrating vectors and α carrying the loadings in each of the r vectors. A preliminary analysis suggest that oil prices have an impact on gas prices with a delay of 13 weeks. Let y t = (g t , b t−τ , T t ) where τ = 91 days. The Johansen test rejects the null of zero cointegrating vectors between the components of y t . The existence of r = 1 cointegrating relation is not rejected and the estimated cointegration vector is, by renormalizing so that the first element be unity,β = (1, −1.0809, 0.0194).
The estimated VECM is as follows. For ease of presentation, unsignificant coefficients, at the 5% level, have been omitted.
It is worthnoting that for the brent prices, no significant linear influence of the past variables is detected. Results not reported here show that the process ( t , ζ t , ξ t ) pass the diagnostic tests for the absence of autocorrelation. Figure 2 displays the series of residuals t for the gas prices. From Figure 3, displaying the empirical autocorrelation function (ACF) and partial autocorrelation function (PACF) of ( t ), it is seen that this series has the characteristics of a white noise. Caution is needed, however, in the interpretation of the significance bands under conditional heteroskedasticity (see Francq and Zakoian, 2009). The same graphs, displayed in Figure 4 for the series ( 2 t ), show that a GARCH effect is present in the data.

Modelling the volatility of gas prices
The volatility models for the series t were estimated over the period april 2000 to december 2004, involving 1192 observations. To have a gauge, the following standard one-regime GARCH(1,1) model was fitted where the standard errors appear in parenthesis. The GARCH coefficients are close to those generally obtained for financial series, with a strong persistence in volatility (α + β = 0.92).
Next, we turn to multi-regimes GARCH(1,1) models, where the regimes are determined by the temperature level. We start by a three-regimes model, where the three classes of temperatures correspond to approximately the same num- ber of observations. This leads to choose s t = 1 when T t < 9, s t = 2 when T t ∈ [9, 14], and s t = 3 when T t > 14, with frequencies in the sample π 1 = 0.35, π 2 = 0.32, π 3 = 0.33. The fitted three-regimes GARCH(1,1) model is as follows.
All coefficients, except the intercept in the first regime, are significant at the 5% level. The most striking point is the difference between the volatility dynamics in the middle regime, compared to the volatilities of the two extreme regimes. The volatility of the second regime is less persistent (α(2) + β(2) = 0.73) with a more convex "news-impact curve". The impact of recent observations on the volatility is stronger than in the low-and high-temperature regimes. It can be noted that the three GARCH(1,1) models are second-order stationary, which entails the global stability with a finite time-dependent variance for t . Note also that the marginal variances within each regimes (ω(j)/(1 − α(j) − β(j)) are roughly the same (around 0.04).
Using the estimated parameters of Model (11) as initial values in the numerical optimization routine, we get the fitted model The effects already noticed for the middle regime (little persistence and strong convexity of the news impact curve) is more pronounced with this five-regimes model. A strong coefficient α is also obtained in the fourth regime. On the contrary, the volatility in all other regimes mainly does not much depend on the last observation. Again, the model is globally stable in the second order sense.
The fitted model is However, many coefficients are found insignificant at the 5% level. Finally, we estimated a model in which the extreme temperatures (low and high) are gathered in the same regime. Letting s t = 1 when T t < 3.2 or T t > 18.5, s t = 2 when T t ∈ [3.2, 9[, s t = 3 when T t ∈ [9, 14[, and s t = 4 when T t ∈ [14, 18.5[, the regimes frequencies deduced from (14) are π 1 = 0.11, π 2 = 0.29, π 3 = 0.32, π 4 = 0.28 (16) and the estimated model is The likelihoods of the different models, displayed in Table 1 allow to compare the different fits. From likelihood ratio tests, at the 5% significance level, • the standard GARCH(1,1) model is not rejected against the 3 regimes model; • the GARCH(1,1) model is however rejected against any model with d > 3; • the 3 regimes model is rejected against the 5 regimes Model (13).
Wald tests not reported here lead to the same conclusions. In the same table, the estimated kurtosis of the variable η t = t /σ t are reported. The biggest kurtosis reduction is obtained with the 5-regimes Model (13).  Table 2 reports Mean-Squared Errors (MSE) of prediction. We re-estimated the different models over the same sample minus the last 500 observations, which were used for the predictions. The estimated models over the sample were very close to those estimated on the whole sample. From the prediction point of view, the 5-regime Model (13) is again the preferred specification.

Conclusion
This paper reviewed a class GARCH models, which allow the volatility to depend on an observed exogenous process. This observability of the state variable makes the model much easier to use than the so-called Markov-switching models in which the regime change is governed by a latent Markov chain. The model can be estimated by QML and a consistent estimator of the asymptotic covariance matrix has been proposed. The methodology has been applied to daily gas prices using the temperature as the exogenous variable. We found evidence of five regimes, with very different dynamics for the volatility in the moderate-temperature regime. The model could be used for prediction purposes, using temperature scenarios. Many extensions, by including more lags in the volatility dynamics or by considering multivariate series, are left for future research. It is hoped that the article will broaden the use of time series models driven by exogenous variables.
where t (θ) is defined as˜ t (θ) withσ t replaced by σ t , we thus havê