1. Introduction
Over the last decade, weather derivatives have emerged as an attractive and interesting new derivative class for risk management and other related fields. Since the purpose of weather derivatives is to allow businesses and other organizations to insure themselves against fluctuations in the weather, a number of insurance, reinsurance companies, banks, hedge funds and energy companies have set up weather trading desks. It is obvious that there are close connections between energy and weather, for example, a natural gas distributor may want to buy weather derivatives to ensure profits in case of a warm winter in which consumers will buy less gas for heating purpose. Hence, weather derivatives will provide a financial instrument for companies or organizations to avoid some unzealous impact of the “bad” weather effects and control the weather risks. Also, additional instruments for hedging in the energy markets require quantitative models of both energy and weather derivatives included.
The weather derivatives market, in which contracts written on weather indices was firstly appeared Over-thecounter (OTC) in July 1996 between Aquila Energy and Consolidated Edison Co. from United States. After that, companies accustomed to trading weather contracts based on electricity and gas prices in order to hedge their price risks realized by weather during the end of 1990s and the beginning of 2000s. Consequently, the market grew rapidly and expanded to other industries and to Europe and Japan. Reported from Weather Risk Management Association (WRMA), an industry body that represents the weather market, recently, the total notional value of the global weather risk market has reached $11.8 billion in last year. With geographic expansion, the OTC market boosted nearly 30% in last year. In this article, we will concentrate on the futures market of temperature derivatives found at the Chicago Mercantile Exchange (CME), which is one of the largest weather derivatives trading platforms. Up to now, the CME has weather futures and options traded based on a range of weather indices for 47 cities from United States, Canada, Europe, Australia and Asia.
As a common sense, weather affects different entities in different ways. In order to hedge these different types of risks, weather derivatives are written on different types of weather variables or weather indices. The most commonly used weather variable is the temperature. Widely used temperature indices include cumulative average temperature (CAT), heating degree days (HDD) and cooling degree days (CDD). They are originated from the energy industry, and designed to correlate well with the local demands for heating or cooling. CAT are defined as the sum of the daily average temperature over the period 
of the contract, the index
where 
is the daily average temperature. It is mainly used in Europe and Canada. In winter, HDD are used to measure the demand for heating, i.e. they are a measure of how cold the weather is and usually used in United States, Europe, Canada and Australia. In contrast, CDD are used in summer to measure the demand of energy used for cooling and a measure of how hot the weather is. They are usually used in United States, Canada and Australia. The definitions for HDD and CDD are given by
and
where the threshold 
takes value of
.
The pricing of a given weather derivative by calculating the expected value of its appropriately discounted payoff is inherently related to weather forecasting and simulation (see Zeng [1] and Cao and Wei [2]). The conventional approach is to identify the probability distribution of the associated index at maturity, and then integrate the payoff of the derivative with respect to it. Dornier and Queruel [3] characterized weather dynamics by mean-reverting Itô diffusions, driven by a standard Brownian motion, and then resorted to Monte Carlo simulation. Davis [4] used the geometric Brownian motion to model the accumulated HDDs (or CDDs). Brody et al. [5] proposed to replace the traditional Brownian motion by a fractional Brownian motion, leading to a fractional Ornstein-Uhlenbeck, allowing incorporate a long memory effect. Benth [6] used the Fourier transformation technique to get the explicit expressions of European and average type of weather contracts, and also showed the pricing formulas satisfy certain Black-Scholes PDE. Benth and Šaltytė-Benth [7] proposed an Ornstein-Uhlenbeck model with seasonal mean and volatility, where the residuals are generated by the class of generalized hyperbolic Lévy process rather than Brownian motion family. The process they used is a flexible class of Lévy process capturing the semi-heavy tails and skewness properties of residuals. A similar model driven by Brownian motion was applied for pricing CAT, HDD and CDD futures in Benth and Šaltytė-Benth [8], they also gave discussions about options written upon these futures. Benth et al. [9] generalized the analysis in the above OU model to higher-order (with lag p) Continuous-time Autoregressive models with seasonal variance for the temperature, and they found the choice of p = 3 turns out to be sufficient to explain the temperature dynamics in Stockholm, Sweden by calibrating the model to more than 40 years of daily observations. More recently, Zapranis and Alexandridis [10] began proposing using a time dependent speed of mean reversion parameter in the Ornstein-Uhlenbeck model and use neural networks to estimate the parameters.
The remaining of this article is organized as follows: In Section 2, we firstly apply existed Lévy driven meanreverting model to the modeling of temperature index of Canadian temperature data, and then present Lévy driven CAR model and its application to Canadian data. In Section 3, we develop the futures pricing techniques with respect to weather futures of CAT, CDD and HDD based on the established two models. Finally, in Section 4, we illustrate the concluding remarks and some possible future research in this direction.
2. Modeling of Weather Index
2.1. Lévy Driven Ornstein-Uhlenbeck Model
As our first proposed model, consider the weather index
, which is the daily average temperature (DAT). We suppose the DAT has a generalization of the dynamics:
(1)
where 
is the seasonal mean level across years given days and 
is the speed in which 
revert to
. 
is assumed to be a measurable and bounded function represents the seasonal volatility of temperature. 
is cádlág, adapted, real valued Lévy process with independent, stationary increments and stochastically continuous. For formal definition of Lévy process, we refer the reader to the introduction paper by Papapantoleon [11]. This model was firstly introduced by Dornier and Queruel [3] with Brownian motion. We will firstly apply this model with Lévy process to Canadian temperature data.
From Equation (1), we can directly use Itô formula for semimartingales to get the explicit solution
(2)
For the purpose of fitting this model to our daily average temperatures, we reformulate Equation (2) by subtracting 
from
,
(3)
with the notation
. The stochastic integral could be approximated by
(4)
Adding 
to both sides
(5)
If we define
, we will have the following three parts for modeling 
step by step:
(6)
where 
is the seasonal component, 
is the cyclical component derived from Equation (5), and 
is the stochastic part.
The data sets we are going to analyse include daily average temperatures (measured in centigrade) from Calgary, AB and Toronto, ON in Canada. The choice of these two cities is intuitively based on the location and importance of these cities. These data are provided by the Canadian National Climate Data and Information Archive (NCDIA) over a period ranging from January 1st 2001 to January 1st 2011, resulting in 3652 records for each city. For simplicity, the measurement on each February 29th was removed from the sample in each leap year to equalize the length of each year.
2.1.1. Linear Trend
Empirically, the daily average temperature should have an increasing trend due to some reasons such as global warming, green-house effect and so on. However, from Figure 1, we cannot verify significant linear trends for both data sets. The slopes for linear fitting of Calgary and Toronto are 
and 
respectively. A potential explanation to this is that our sample size is reasonably small with a period of 10 years data. Larger sample size will show us relatively more significant increasing trend using least square fit. Hence, at this stage we will ignore the linear trend, and keep a constant mean level in the seasonal trend.
2.1.2. Seasonal Trend
In order to model the seasonality pattern, we use the trigonometric function with the form
(7)
where 
is the mean level, 
is amplitude of the mean. The period of the oscillations is one year (365 days neglecting leap year), so the weight
. 
represents the phase angle. The results of parameters estimation using least square method are reported in Table 1.
2.1.3. Cyclical Component
To remove the adjacent correlation in the time series, make the time series stationary, we model the cyclical component by recursive regression based on Equation (5). The cyclical component
can be modeled by:
(8)
The results of
estimated implicitly by an autoregressive model and corresponding
are given in Table 2.
The autocorrelation function (ACF) plots for the squared residuals in Figure 2 clearly reveal seasonality of correlation for both cities. In the squared residuals ACF plots, the autocorrelation violated 
interval for lots of lags, and the wave for the ACF is a sign of
Figure 1. Linear trends of Calgary and Toronto.
Table 1. Estimated parameters for seasonal function.
Table 2. Estimated parameters in regression.
seasonal pattern.
2.1.4. Seasonal Volatility
We consider the seasonal pattern in the residual after removing the cyclical component as the seasonal volatility. The motivation behind seasonal volatility is straight forward. The temperature appears bigger variation in the winter and smaller variation in the summer. By Equation (5), the random term has the form
. To get the seasonal volatility
, we have the following methods:
• Organize the data into 365 groups, one for each day of the year;
• Find the means of the squared data in each group. These means are crude estimation for the expected squared data.
;
• To smooth the estimation, first way is trying the smoothing technique. In details, we could take logarithm of the estimation and use moving average method to smooth the crude estimation. Then the smoothed estimation of 
could be get by ex-
Figure 2. Empirical ACF after regression.
ponentiation of the result. Alternatively, use the truncated Fourier series (9) to fit the estimation from previous two steps. The Fourier series will indicate the presence of different cycles in 
as a deterministic function.
(9)
Figure 3 shows the fitted and smoothed 
by three days moving average technique. More average smoothing range will make the estimation smoother, but for our Canadian temperature study, the variation of the temperature in more days would be relatively larger than that of three days. Hence, the choice of three days moving average makes sense for both the smoothing and accuracy purpose.
After removing the estimated seasonal volatility, the ACF for the data and squared data are plotted in Figure 4. It shows that the seasonality pattern is completely removed from the squared data.
On the other hand, we also tried the way of fitting the crude seasonal volatility by truncated Fourier series with other three. The values of the estimated parameters are given in Table 3.
Corresponding fitted results are shown in Figure 5.
Again, the seasonal pattern in the ACF plot was removed from Figure 6. So, we could find that both methods can model the seasonal volatility appropriately.
2.1.5. Random Noise
Then we head to fit the random part 
remained from
Table 3. Truncated fourier series parameters.
above. The classical literatures of temperature modeling suggest the random noise are independent and identically normal distributed, (see Dornier and Queruel [3], Alaton et al. [12]). However, we could see from the QQ-plot and normal fit of the random noise for our data, the normal distributed random noise may lead misdescription to the data. Moreover, modeling of random noise using Gaussian process is not always the case.
Instead, we try to apply the generalized hyperbolic family process to model the random part for our data, based on Benth and Šaltytė-Benth [7]. As the name suggests, it is super-class of processes such as normal-inverse Gaussian, variance-gamma, hyperbolic and so on. This will allow the family better capture the features of the data set and relatively easy to analyze from the Lévy properties it holds.
The generalized hyperbolic distribution is an infinitely divisible distribution which can be defined with the density function
Figure 3. Empirical and smoothed seasonal variation.
Figure 4. ACF after removing seasonal variation.
where 
is modified Bessel function of the third kind with index
. Parameters of the generalized hyperbolic distribution can be estimated by maximum likelihood method. In detail, for a vector of observations
, the maximum likelihood estimate of the parameters 
can be obtained by maxi-
Figure 5. Fit of volatility by truncated Fourier series.
Figure 6. ACF after removing seasonal variation.
mizing the log-likelihood function for generalized hyperbolic distribution:
McNeil et al. [13] described a modified EM method, which is called multi-cycle, expectation, conditional estimation (MCECM) algorithm to optimize the log-likelihood function. We use this algorithm written as “ghyp” package in “R” to estimate parameters. Estimated parameters are presented in Table 4.
In order to check the goodness of fitting, we simulated a series of generalized hyperbolic distributed random
Table 4. Estimations of generalized hyperbolic parameters.
variables using estimated parameters, and plotted the empirical cumulative distribution function, then compared with our residual observations’ empirical cumulative distribution function plot. The results are shown in Figure 7. The result implies our calibration to the model is reasonably fair for Canadian data.
2.2. Lévy Driven Continuous-Time AR Model
Another suitable class of stochastic process we could use to analyze the evolution of temperature is an extension of the previous OU process to the multi-factor OU process. The model is called Continuous-time Autoregressive (CAR) model, which is also a subclass of the general Continuous-time Autoregressive Moving-average (CAR-MA) model introduced by Brockwell and Marquardt [14]. The intuition behind using this model is the temperature memory is consistent with high-order AR model, and the seasonality of mean and volatility is also involved in the model. In other words, the CAR model combine the mean-reverting property of OU process and the memory of AR model together to better capture the features of temperature.
Followed by the classical stochastic assumption, suppose given a probability space
. Let 
be a stationary solution to the stochastic differential equation in 
for 
(10)
where the operator
denotes the differential operator with respect to
, 
is the p-order differential operator and 
is 1-dim Lévy process. Let 
be a column vector defined as
We could also rewrite the Equation (10) in the statespace form.
(11)
(12)
The vector 
is usually recognized as state vector. The standard deviation of the noise is described by a continuous function
, we usually refer the
Figure 7. CDF of simulated series and observations.
function 
as the volatility of the process
. The constant 
denotes the order of the CAR model.
Assume additionally the seasonal function 
is bounded and continuously differentiable on
. The temperature dynamics is formulated by:
(13)
where 
is the first entry of the state vector 
in Equation (11). Note that for Equation (12), since we have the function 
to capture the deterministic temperature level, the OU process should ideally revert to zero. This is the main reason we have only the meanreverting rate 
but without level in dynamic of 
Again, by using the Itô lemma for semimartingale, the explicit solution for the stochastic differential Equation (12) has the form
(14)
If we construct the expectation of the Wiener-Lévy integral in Equation (14) equals to zero (it is obviously true for the Itô integral in the special case), once all eigenvalues of the matrix 
have negative real parts, the process 
then have a mean equal to zero. We could still use simple intuition to have the idea, when eigenvalues of the matrix 
have negative real parts and take 
then the expectation goes to zero. This result implies the temperature dynamic on average will mainly depend on the seasonal function 
by Equation (13).
Notice that we could have a link point of view to this CAR model from the previous OU model. Consider CAR(1) model as an example. In this case, 
in Equation (12), the matrix 
will simply becomes a constant
. The process 
and by Equation (14), the temperature dynamics
Our proposed CAR model is therefore consistent with the single-factor OU model and also satisfies the requirements of regression and mean-reverting properties. For simplicity of estimation, we connect our CAR model to the classical AR model. Benth et al. [15] (Lemma 10.2) provided an explicit way to find the relation between the coefficients 
in CAR model and 
in AR model by induction as
(15)
with coefficients 
and initially 
Next, we will recap the calibration procedure to the same empirical Canadian daily average temperature data based on the CAR model.
2.2.1. Linear and Seasonal Trend
As we did for OU model in Section 2.1.2., we could combine the linear and seasonal trend in the OU model to construct the seasonal function
in CAR model as
(16)
The coefficients for trigonometric function in Equation (16) play the same role in Equation (7). Followed by the description of linear trend in Section 2.1.2., we eliminated the linear part in 
again for CAR model, and the fitted results for seasonal part are the same to that in Table 1.
2.2.2. Cyclical Component
The PACF plots in Figure 8 clearly show that the evolution of detrended temperature 
could be modeled by the classical AR(3) model. Motivated by the PACF plots, by Equation (15), we have 
The estimated parameters are outlined in Table 5.
PACF plots in Figure 9 show the significant correlation in the first three lags is removed by regression, but still a distinct seasonality in the ACF plots of squared residuals.
2.2.3. Seasonal Volatility
For the CAR model, we will apply the truncated Fourier series in Equation (9) to model the seasonal volatility. We choose the degree N = 1 in Equation (9), since there is no significant improvement for the result by adding degree. The fitted parameters are shown in Table 6.
In Figure 10 we present the empirical 
and fitted
, both of them show that the fluctuation in the cold season are considerably larger than that of the mild season.
Figure 11 tells us that seasonal pattern in the squared residual is removed successfully.
Table 5. Estimations of cyclical component.
Table 6. Truncated Fourier series parameters.
2.2.4. Random Noise
We try again apply the generalized hyperbolic process to model our random noise part. The method is the same to that for OU model. Table 7 shows the fitted parameters for our Canadian cities data.
Also as we did for OU model, we simulated a series of generalized hyperbolic random variables with estimated parameters, then compare the empirical CDF with that of historical data in Figure 12. The empirical CDF coincides with simulated CDF almost perfectly, so we may expect this model fairly well describe our Canadian temperature data.
Figure 10. Empirical and fitted seasonal volatility.
Figure 11. Empirical ACF plot of squared residuals.
3. Pricing of Weather Derivatives
In this section, we will derive the temperature futures prices written on CAT, CDD and HDD, which constitute the three main classes of future products at CME market.
3.1. Future Pricing of Lévy Driven OU Model
Consider the price dynamic of future written on CAT over specific time period 
with 
Firstlyassume the daily average temperature follows stochastic differential Equation (1) with 
being Lévy process and constant continuously compounding interest rate 
The future price at time 
based on CAT under risk-neutral probability measure 
will be 
adapted stochastic process 
satisfying
(17)
Figure 12. Empirical ACF plot of squared residuals.
Table 7. Estimations of generalized hyperbolic parameters.
By adaption of the stochastic process 
with respect to measure 
the future price
(18)
To derive the explicit future price in Equation (18), we need to specify a risk-neutral probability measure 
However, the commodity market is typical incomplete market, since most of commodity trades impose big transaction and storage cost. For our case, the underlying temperature is even not possible to be stored and traded. These features break down the classical hedging approach used to derive the unique fair price of derivatives. Because of the incompleteness of the temperature market, any probability measure Q being equivalent to the real measure P is a risk-neutral measure. Next, following the analysis in Benth and Šaltytė-Benth [16], we will specify a class of risk-neutral measure via Esscher transform. Before using the Esscher transform, we need an integrable condition for the Lévy measure to ensure the existence of moments of underlying asset process.
Condition 1. [Benth et al. ([15], p. 74)] 
such that the Lévy measure
satisfies the integrable condition
(19)
almost surely.
The constant 
determines the order of the moment that is finite for the underlying process. With this condition, we recap the following results from several sources that we will use in our further analysis (see Section 3.2).
Lemma 1. [Benth and Šaltytė-Benth [16]] Let function 
be a bounded and measurable function and Condition 1 holds for 
then
(20)
where 
is the cumulant function of the Lévy process 
is the Lévy symbol defined by map 
such that the characteristic function of the Lévy process 
Now assume 
is a measurable and bounded function. Let us consider the stochastic process
(21)
The probability measure 
can be defined by the Esscher transform
(22)
where 
is the indicator function over a probability space
. Obviously this measure is equivalent to the real measure 
and transformed from the probability measure 
so an equivalent risk-neutral measure. A way to get a more flexible class of equivalent martingale measure is to use time varying 
when fitting the observed forward curve. Note that the market price of risk is represented by 
in this case.
Then, consider the existence of the fair forward price dynamics 
given by Equation (18). We need to have another condition such that
is finite almost surely.
Condition 2. [Benth and Šaltytė-Benth [16]] Assume a constant 
satisfies 
such that Condition 1 holds, then the price dynamics
almost surely.
With the existence of explicit risk-neutral probability measure and fair pricing technique, we could derive the following results for the future price written on CAT.
Lemma 2. [Benth and Šaltytė-Benth [7]] Given a measurable and bounded function 
the expectation of Lévy integral under risk-neutral measure over 
is
With these lemmas, we could get the explicit formula for the future price written on CAT.
Theorem 1. [Benth and Šaltytė-Benth [7]] The future price 
at time 
written on CAT index is given by
Next, for pricing of CDD and HDD futures, recall that the CDD and HDD over a measurable time 
is
Before introduce the pricing technique for the CDD and HDD future, we should notice that there is a parity relation between futures prices of CAT, CDD and HDD. This relation will help us to get one of future price of CDD or HDD, once we have the CAT future price and the other one of them.
Theorem 2. [Benth et al. ([15], p. 279)] The relation between futures prices of CAT, CDD and HDD is
(23)
With this theorem, we will focus on the pricing of CDD future, but the pricing technique for CDD future is also analogously applicable to HDD future.
Based on risk-neutral pricing theory, the CDD future price satisfies
with risk-free interest rate 
and risk-neutral measure 
Given the risk-neutral measure derived from Esscher transform, since the CDD future under the measure is adapted, the risk-neutral CDD future price is derived as
(24)
Theorem 3. [Benth and Šaltytė-Benth [7]] For
the logarithm of characteristic function 
of 
under measure 
is given by
(25)
where 
is the cumulant generating function of 
This theorem provides us a way to get the characteristic function of Lévy integral of the type 
under the risk-neutral measure 
With the characteristic function, the pricing approach we could try is using the celebrated Fourier inversion theorem to get the numerical probability density function of 
under 
measure, and then using the density function to price CDD future.
Theorem 4. (Fourier Inversion Theorem) [Hewitt and Stromberg ([17], p. 409)] Let 
denote the cumulative density function of random variable 
the corresponding probability density function 
is integrable in Lebesgue sense, i.e. 
The characteristic function of 
is defined as
Then,
(26)
To calculate the probability density function from characteristic function, consider the density 
and the discrete Fourier transform (DFT)
We could use the trapezoid rule to approximate the integral in Equation (26).
(27)
where 
The last step of approximation is because when we take 
to be large enough, 
The purpose to have the upper bound of summation as 
is to match with the form of DFT.
Now we propose an approach to calculate the probability density function analytically based on the construction of Chourdakis [18]. Note that the application of DFT will result in a set of integral approximation which have the form of Equation (27) based on
Let 
then 
For the return values 
these values are setting to be also equidistant with the grid spacing 
By contrasting Equation (27) and DFT, we have 
Then set the return grids
(28)
where 
is a parameter to control the return range. For the center returns around zero, we could set 
With these settlements, from Equation (27), the probability density function has the form
(29)
where 
This summation has the form of DFT, could be computed efficiently.
Theorem 5. Numerically, the CDD future
where 
is given by Equation (29) and 
are constructed by Equation (28) in the approach of fast Fourier transform.
Proof. Recall the CDD future is priced by Equation (24). By Fubini theorem, the definition of expectation and property of Riemann integral,
3.2. Future Pricing of Lévy Driven CAR Model
Now we assume the temperature dynamic follows Lévy driven CAR model, introduced by Equations (11) and (12). The risk-neutral measure 
is given by Equation (22), and the Existences of futures prices for CAT, CDD and HDD are preserved under Condition 2.
Based on the derivation of futures prices of CAT, CDD and HDD indices under Lévy driven OU model, we have the following results for the futures prices of these indices under the Lévy driven CAR model.
Lemma 3. The cumulated temperature on the time interval 
under Lévy driven CAR model is
Proof. Integrating the original Equation (12) over the time interval 
we have
Then, by using the solution of Lévy stochastic model, Equation (14), Equations (11) and (13),
So, the result is proved. 
Theorem 6. The future price of CAT
Proof. By Equation (18), we could use Lemma 3 to expand each term separately, and then use Lemma 2 to get the final explicit formula. 
For CDD or HDD futures pricing, by using Theorem 3the characteristic function of 
could be found as 
where
Then, followed by the approach in the OU model, we have Theorem 7. Numerically, the CDD future price
where 
is given by Equation (29) and 
are constructed by Equation (28) in the approach of fast Fourier transform to compute 
The proof of this theorem is similar to that of Theorem 5.
4. Summary, Conclusion and Future Work
This paper investigates the problem of modeling and pricing weather derivatives and their application to Canadian data. We have shown that two models, including Ornstein-Uhlenbeck model and Continuous-time Autoregressive model, both driven by Lévy process are fairly suitable to capture the evolution of Canadian cities temperature. Also, based on these two models, we derived approaches for risk-neutral prices of future contracts written on CAT, CDD and HDD.
For the Ornstein-Uhlenbeck model, we followed the analysis developed in Benth and Šaltytė-Benth [7], and calibrated the model to empirical study of our Canadian temperature data. As an extension, we derived a numerically explicit pricing formula for futures written on CDD and HDD futures. For the Continuous-time Autoregressive model, we extended the original model introduced by Benth et al. [15] to the one with Lévy process as the random variable instead of Brownian motion. And we also, calculated the explicit formulas for futures written on CAT, CDD and HDD indices under the structure. We would like to mention that our approaches and results could be applied to other temperature data set rather than Canadian ones with some modifications and adjustments.
This work can be extended in several ways. For instance, based on our empirical studies of Canadian data, the volatility possesses significant seasonal behaviour under the existed models. One could consider using regime-switching techniques to model this seasonal volatility under the structure of the established models in this paper. For instance, depending on the local appearance of the weather in Canadian cites, one can intend to use twostate Markov chain for modeling of Calgary’s volatility and four-state one for Toronto’s.
As an application to risk management in industries, we could also consider the dynamic hedging strategy of other futures by using weather future. Take the energy futures for example, if we consider a portfolio at time 
containing one unit of energy (e.g. heating oil) future 
and 
units of weather futures 
both with maturity at time T. Assume the portfolio has value 
at time 
then
If the portfolio is self-financing, the change in this portfolio in a small amount of time 
is given by
Assume the energy future and weather future followed by a model we have built in this paper, the future hedge ratio
In this case, the stochastic component of portfolio vanishes and the portfolio value is hedged. Thus, in order to hedge an energy future, we can short 
shares of weather futures and the portfolio is hedged dynamically.
Future discussion about this idea will be centered on two viewpoints. Firstly, further research need to be done for the dynamics of energy future and weather future. As such, explicit dynamic hedging conclusion regarding to hedging risks of energy commodities could be founded. The other point of view is allowing other closely related derivatives to be involved into the constructing of portfolio. For instance, we could also take the newly burgeoning carbon dioxide derivatives into account in the portfolio. These considerations will be a comprehensive tool to integrate into hedging of energy commodities’ risks.