There are many management methods and problems about energy risk. This paper focuses on how to use energy futures to hedge the risk of energy spots. It is assumed that the rates of return on energy futures and spots are r 1 and r 2 , respectively; the mean and standard deviation of r 1 are μ 1 , σ 1 ; the mean and standard deviation of r 2 are μ 2 , σ 2 . The portfolio rate of return r p formed by 1 unit of long energy spots and h units of short energy futures can be expressed as: r p = r 1 − h r 2 (1)

The mean and variance of r p are μ p = μ 1 − h μ 2 and σ p 2 = σ 1 2 + h 2 σ 2 2 − 2 h ρ σ 1 σ 2 respectively, in which ρ is the correlation coefficient of rates of return on energy spots and energy futures. Risk hedging strategy focuses on finding the best hedge ratio h to minimize the risk of the hedge portfolio. According to the different risk measurement indicators, risk hedging strategies of futures and spots are different. The classic and most commonly-used risk measurement is the variance proposed by Markowitz (1952). The minimum variance hedging strategy is to find the hedging ratio that minimizes the variance of the portfolio.

It is assumed the cumulative distribution function of the rate of return r p on hedging portfolio is F ( x , h ) , and at a given confidence level of 1 − α , the value-at-risk V a R ( h ) of the hedging portfolio is expressed as (Jorion, 2000): V a R ( h ) = − inf { x ∈ R : F ( x , h ) ≥ α } (2)

The confidence level 1 − α is given in advance, usually by investors according to their own preferences or by regulatory agencies. In the condition that the distribution function F ( x , h ) satisfies the continuity, V a R ( h ) is the opposite of the lower α quantile of r p . Although VaR is a widely-used risk measurement index, many scholars show the defect that VaR does not satisfy the subadditivity and it ignores the tail risk (Acerbi and Tasche, 2002). Therefore, Rockafellar and Uryasev (2000) proposed the concept of CVaR, which refers to the mathematical expectation of all losses exceeding the VaR level. The mathematical expression of the CVaR of the hedging portfolio is: C V a R ( h ) = − E [ r p | r p ≤ − V a R ( h ) ] (3)

Risk hedging strategy on the basis of CVaR focuses on finding the hedging rate h that minimizes the C V a R ( h ) of the hedge portfolio. Under the assumption that r p is subject to normal distribution, the corresponding C V a R ( h ) can be expressed as the linear function of the mean value and standard deviation (Alexander and Baptista, 2004), namely: C V a R ( h ) = φ ( z α ) α σ p − μ p (4) where σ p is the standard deviation of the combined rate of return. z α is the lower α quantile of the standard normal distribution. φ is the density function of the standard normal distribution, The normal distribution assumption is not true in the real market. An equivalent definition of CVaR given by Rockafellar and Uryasev (2000) without any distribution setting is: C V a R ( h ) = m i n v ∈ R F α ( h , v ) (5) where, F α ( h , v ) = v + α − 1 E [ ( − r p − v ) + ] , and ( x ) + = m a x ( x , 0 ) . Meanwhile, the authors pointed out that the minimum CVaR model is equivalent to the following optimization problem: m i n h ∈ R C V a R ( h ) = m i n ( h , v ) ∈ R × R F α ( h , v ) (6)

Given the distribution function or density function of r p , F α ( h , v ) can be obtained. However, we seldom know the specific form of distribution function or density function in real life and need to estimate it by using sample data. Supposing { r 1 , t } t = 1 T and { r 2 , t } t = 1 T are sample data of return rate on spots and futures, respectively, and { r p , t } t = 1 T is the sample data of portfolio rate of return in the period of T , and r p , t = r 1 , t − h r 2 , t , thus the kernel estimator of the density function of the combined rate of return is: f ^ ( x , h ) = 1 T b ∑ t = 1 T g ( x − r p , t b ) (7) where G ( z ) = ∫ − ∞ z g ( u ) d u , and g ( z ) is a kernel function selected by the researcher. According to the study by Li and Racine (2007), the Gauss kernel function is a good candidate if we intend to estimate the density function and distribution function of the univariate, that is, g ( z ) = ( 2 π ) − 1 ⁡ e x p ( − z 2 / 2 ) . b is bandwidth. According to the study by Li and Racine (2007), we set b = 1.06 × T − 1 / 5 × σ ^ p = a σ ^ p , in which a = 1.06 × T − 1 / 5 , and σ ^ p is the sample covariance of the portfolio rate of return. Thus, under the kernel estimation framework, according to the kernel estimator (7) of density function, the kernel estimator of F α ( h , v ) can be expressed as: F ^ α ( h , v ) = v + α − 1 E ^ [ ( − r p − v ) + ] = v + α − 1 ∫ − ∞ + ∞ ( − x − v ) + f ^ ( x , h ) d x = v + α − 1 ∫ − ∞ + ∞ ( − x − v ) + 1 T b ∑ t = 1 T g ( x − r p , t b ) d x = v + 1 T b α ∑ t = 1 T ∫ − ∞ − v ( − x − v ) g ( x − r p , t b ) d x = v + 1 T α ∑ t = 1 T ∫ − ∞ R t ( − b y − r p , t − v ) g ( y ) d y = v − 1 T α ∑ t = 1 T ( ( r p , t + v ) G ( R t ) + b H ( R t ) ) (8) where y = ( x − r p , t ) / b , R t = ( − v − r p , t ) / b and H ( x ) = ∫ − ∞ x y g ( y ) d y . The minimum CVaR model under the kernel estimation framework can be expressed as: m i n h ∈ R C V a R ^ ( h ) = m i n ( h , v ) ∈ R × R F ^ α ( h , v ) = m i n ( h , v ) ∈ R × R v − ( T α ) − 1 ∑ t = 1 T ( ( r p , t + v ) G ( R t ) + b H ( R t ) ) (9)

Similar to the study of Yao et al. (2013), the convexity theorem of optimization problem (9) is given as follows:

Theorem 1

Optimization Problem (9) is a convex optimization problem.

Proof: Because the constraint set of the optimization problem (9) is a nonempty convex set, we only need to prove that the Hessian matrix of the objective function F ^ α ( h , v ) is a positive semi-definite matrix.

According to the function F ^ α ( h , v ) and the formula R t = ( − v − r p , t ) / b , we have:               ∂ F ^ α ( h , v ) ∂ R t = − 1 T α ∑ t = 1 T ( ( r p , t + v ) g ( R t ) + b R t g ( R t ) ) = 0 (10) ∂ R t ∂ v = − 1 b , ∂ R t ∂ h = r 2 , t b + v + r p , t b 2 ∂ b ∂ h (11) Furthermore, we take first order derivatives of v and h by the function F ^ α ( h , v ) . By simplifying the formula using the results of (10) and (11), we can obtain: ∂ F ^ α ( h , v ) ∂ v = 1 − 1 T α ∑ t = 1 T G ( R t ) (12) ∂ F ^ α ( h , v ) ∂ h = 1 T α ∑ t = 1 T [ G ( R t ) r 2 , t − H ( R t ) ∂ b ∂ h ] (13)

Furthermore, by Formulas 12, 13, we have the second partial derivatives, as follows: ∂ 2 F ^ α ( h , v ) ∂ v 2 = 1 T b α ∑ t = 1 T g ( R t ) (14) ∂ 2 F ^ α ( h , v ) ∂ v ∂ h = − 1 T α ∑ t = 1 T g ( R t ) ∂ R t ∂ h (15) ∂ 2 F ^ α ( h , v ) ∂ h ∂ v = 1 T α ∑ t = 1 T [ g ( R t ) ∂ R t ∂ v r 2 , t − R t g ( R t ) ∂ R t ∂ v ∂ b ∂ h ] = − 1 T α ∑ t = 1 T g ( R t ) [ r 2 , t b + v + r p , t b 2 ∂ b ∂ h ] = − 1 T α ∑ t = 1 T g ( R t ) ∂ R t ∂ h (16) ∂ 2 F ^ α ( h , v ) ∂ h 2 = 1 T α ∑ t = 1 T ( g ( R t ) ∂ R t ∂ h r 2 , t − R t g ( R t ) ∂ R t ∂ h ∂ b ∂ h − H ( R t ) ∂ 2 b ∂ 2 h ) = 1 T α ∑ t = 1 T ( g ( R t ) ∂ R t ∂ h ( r 2 , t − R t ∂ b ∂ h ) ) − 1 T α ∑ t = 1 T ( H ( R t ) ∂ 2 b ∂ 2 h ) = 1 T α ∑ t = 1 T ( g ( R t ) ( ∂ R t ∂ h ) 2 b ) − 1 T α ∑ t = 1 T ( H ( R t ) ∂ 2 b ∂ 2 h ) (17)

Then, the Hessian matrix of the objective function F ^ α ( h , v ) can be expressed as: Θ = ( ∂ 2 F ^ α ( h , v ) ∂ v 2 ∂ 2 F ^ α ( h , v ) ∂ v ∂ h ∂ 2 F ^ α ( h , v ) ∂ h ∂ v ∂ 2 F ^ α ( h , v ) ∂ h 2 ) = 1 T α ( ∑ t = 1 T g ( R t ) 1 b − ∑ t = 1 T g ( R t ) ∂ R t ∂ h − ∑ t = 1 T g ( R t ) ∂ R t ∂ h ∑ t = 1 T ( g ( R t ) ( ∂ R t ∂ h ) 2 b ) )           − 1 T α ( 0 0 0 ∑ t = 1 T ( H ( R t ) ∂ 2 b ∂ 2 h ) ) (18) Θ 1 = 1 T α ( ∑ t = 1 T g ( R t ) 1 b − ∑ t = 1 T g ( R t ) ∂ R t ∂ h − ∑ t = 1 T g ( R t ) ∂ R t ∂ h ∑ t = 1 T ( g ( R t ) ( ∂ R t ∂ h ) 2 b ) ) = 1 T α ∑ t = 1 T ( g ( R t ) ( 1 b − ∂ R t ∂ h b ) ( 1 b − ∂ R t ∂ h b ) ) (19)

Kernel function g ( ⋅ ) > 0 , bandwidth b > 0 , and 1 T α > 0 . Therefore, Θ 1 , the sum of T positive semi-definite matrixes, is still a positive semi-definite matrix. Θ 2 = − 1 T α ( 0 0 0 ∑ t = 1 T ( H ( R t ) ∂ 2 b ∂ 2 h ) ) (20)

For any z ∈ R , H ( z ) = ∫ − ∞ z y g ( y ) d y , if z < 0 then y g ( y ) < 0 for any y ∈ ( − ∞ , z ) , since g ( y ) > 0 . If z ≥ 0 , we have H ( z ) = ∫ − ∞ z y g ( y ) d y = ∫ − ∞ − z y g ( y ) d y + ∫ − z z y g ( y ) d y     = ∫ − ∞ − z y g ( y ) d y < 0 (21)

Because y g ( y ) is an odd function. Therefore H ( z ) ≤ 0 for any z ∈ R .

Because − 1 T α < 0 , ∂ 2 b ∂ 2 h = a σ ^ 1 2 σ ^ 2 2 σ ^ p 3 ( 1 − ρ ^ 2 ) ≥ 0 and H ( z ) ≤ 0 , so we prove Θ 2 is a semi-positive definition matrix immediaely. Thus, optimization Problem (9) is a convex optimization problem. There are many algorithms to solve the convex optimization problem. In particular, for the one-dimensional optimization problem, we can obtain its global optimal solution easily.

First, by using the nonparametric method, we measure the risks of the two products in the whole samples, with the result noted as NPK-CVaR. Meanwhile, the measurement result of CVaR under normal distribution is given as Norm-CVaR for comparison. In addition, because values of CVaR are related to the level of confidence, a series of different confidence levels are taken to test the measurement results. The results are shown in Table 2.

Table 2 shows that with all the four products and confidence levels 1 − α , Norm-CVaR is less than NPK-CVaR, indicating that the normal distribution assumption always underestimates the risk of energy spots and futures in reality. In addition, according to the comparison between spots risk and futures risk, the measurement results of crude oil and gasoline show that no matter under nonparametric estimation method or normal distribution assumption, the risk value of futures are greater than that of the corresponding spots.

To see the differences of risks measured by the two different methods more directly, we use kernel estimation method and normal distribution to fit the sample distribution of the rates of return on the four products (Figure 1). Specifically, based on the sample data of the rates of return on the four products, the density function of the real data is estimated through kernel estimation, and the image is noted as NPK. As a comparison, the normally distributed density function is used to fit the density function of the sample data. That is to say, supposing that the real data are in normal distribution, we estimate the mean and the standard deviation through sample data, and the normally distributed density function can be obtained and noted as Norm in the figure. It can be seen from the four figures that compared with the density function in normal distribution, density functions estimated by NPK have higher peaks and thicker left tail and right tail. This characteristic of “leptokurtosis and fat-tail” is vital for measurement and management of energy risk. For the long position (short position), the thicker left tail (right tail) indicates a greater probability of larger losses in the real energy market, which is seriously underestimated under normal distribution. This is consistent with the conclusion of Table 2; that is, the risk measurement under the normal distribution assumption tends to underestimate the risk.

This section uses the futures data and spots data of crude oil and gasoline to study the risk hedging based on CVaR. Specially, we use the futures of the crude oil and gasoline to hedge the risk of the corresponding spots. As a comparison, we present the hedging performance under the NPK CVaR and the Norm CVaR. The previous analysis has shown that the Norm CVaR tends to underestimate the actual risk, and CVaR can be expressed as a linear function of the mean and standard deviation under normal distribution assumption, while only taking into account the information on the first two order moments. When the financial market data are in non-normal distribution, investors will pay more attention to higher moments of market data in addition to the first two moments, such as mean and variance. Therefore, Cao et al. (2010) proposed a semi-parametric method to estimate CVaR. In fact, their method is the Cornish-Fisher expansion of CVaR (CFM CVaR), which includes the information on the third moment and fourth moment of the market data without any distribution assumption. This method is an improvement of normal method. Although CFM CVaR reflects the information on the first four moments of the risk factors, the information on higher order moments cannot be included because of the complexity of the formulas. The nonparametric kernel estimation method can get the sample distribution function of the market data exactly without any distribution setting, and therefore that the information on higher moments can be obtained. Theoretically, if the sample size is large enough, the nonparametric kernel estimation method can get the information on all moments of the risk factors. Therefore, we compare the performance of the three methods in the actual risk hedging. We consider both the static hedging strategy and the dynamic hedging strategy. Static hedging strategy, which is the optimal hedging strategy obtained from the training subsamples, is applied to the training subsamples and testing subsamples to test the in sample and out of sample performances of the hedging strategies.

It is a defect that the static hedging assumes that the distribution of the sample data in the sample interval is unchanged, while the actual financial market environment is constantly changing. Therefore, if the sample interval is too long and the distribution of rate of return changes, then there may be deviation in the static hedging. In dynamic hedging, long intervals are divided into several short intervals, in which the distribution of rates of return is assumed to be unchanged. This avoids the problem in static hedging that the distribution of the rates of return is unchanged in a long interval. However, the defect in dynamic hedging is that the position of the futures needs to be constantly adjusted, which increases the cost of risk hedging. The dynamic hedging strategy is constructed as follows: the optimal hedge ratio is obtained from the initial 1565 training subsamples, and is applied to the next 10 trading days (2 weeks). The training subsamples are then updated by replacing the oldest 10 samples by the data of the 10 trading days, as mentioned earlier, while the training sample size of 1565 is unchanged. The new hedge ratio is obtained from the updated training subsamples, and is applied to the next 10 trading days, and so on until to the end of the sample period. This method reflects the latest transaction information by dynamically adjusting the training subsamples. Finally, based on the dynamic hedge ratio, the sequence of the portfolio rate of return is obtained and the out of sample performance (e.g., mean, standard deviation and CVaR) are calculated, see Table 3.

Table 3 shows the mean, the decline in the standard deviation (%ΔStd), and the decline in the CVaR (%ΔCVaR) of the portfolio rate of return after hedging based on the three methods. After being hedged by three methods, the CVaR and the standard deviation declined greatly. This shows whether the decline in the CVaR is the largest based on NPK CVaR from the point of view of a single asset or an average value. Specifically, NPK CVaR is superior to Norm CVaR, and Norm CVaR is superior to CFM CVaR. This indicates that as a modification to Norm, semi-parametric CFM is not satisfactory in practice. The decline in the standard deviation is the greatest in the Norm CVaR method, mainly because CVaR is the linear function of the standard deviation in the Norm method, and to minimize Norm CVaR is to minimize the standard deviation. The standard deviation is a symmetric index. Its decline may be caused by the decline, either in the left tail or in the right tail of the distribution. This shows that the decline in the standard deviation in Norm method is larger than that in NPK method, which is largely due to the decline in the right tail in the table because the left decline of CVaR in Norm is smaller than that in NPK. Similarly, the sequence of standard deviation declines from large to small is in the sequence of Norm, NPK, and CFM. Finally, the conclusion by the mean is mixed. The means in the three methods after hedging are all positive in sample, while they are negative out of sample. This happens because the mean of the rate of return on spots before hedging is positive in sample, while it is negative out of sample.

From Figures 2–5, the histograms of the pre-hedging and post-hedging rates of return on the two energies in and out of samples are shown, respectively. Only the hedging image in NPK CVaR is given, because the risk decline in it is the largest. It can be seen directly from the figures that after hedging, the histograms are more concentrated, and the tail samples at both ends are less with most of the samples distributed on the right side of zero. This shows that hedging significantly reduces the risk in the left tail. The empirical results show that the NPK method that we proposed is more effective to measure the actual energy risk and carry out more effective risk hedging.