Modelling and Computation in the Valuation of Carbon Derivatives with Stochastic Convenience Yields

The anthropogenic greenhouse gas (GHG) emission has risen dramatically during the last few decades, which mainstream researchers believe to be the main cause of climate change, especially the global warming. The mechanism of market-based carbon emission trading is regarded as a policy instrument to deal with global climate change. Although several empirical researches about the carbon allowance and its derivatives price have been made, theoretical results seem to be sparse. In this paper, we theoretically develop a mathematical model to price the CO2 emission allowance derivatives with stochastic convenience yields by the principle of absence of arbitrage opportunities. In the case of American options, we formulate the pricing problem to a linear parabolic variational inequality (VI) in two spatial dimensions and develop a power penalty method to solve it. Then, a fitted finite volume method is designed to solve the nonlinear partial differential equation (PDE) resulting from the power penalty method and governing the futures, European and American option valuation. Moreover, some numerical results are performed to illustrate the efficiency and usefulness of this method. We find that the stochastic convenience yield does effect the valuation of carbon emission derivatives. In addition, some sensitivity analyses are also made to examine the effects of some parameters on the valuation results.


Introduction
As we all have known, climate change is one of the main environmental problems. Over the last decades, numerous scientific studies have proved that greenhouse gases such as carbon dioxide undoubtedly contribute to the climate change a lot. The accumulation of greenhouse gases, particularly CO 2 , contributes to the high concentration of solar energy in the air, and the greenhouse effects can be reflected by the increase of the number of extreme weather events, such as tsunamis, floods and droughts. The greenhouse gas emission increases so fast that the temperature in the atmosphere is rising, which is now a big threat to all the species on the earth, and leads to the glacier melting and the sea level rising. Thus, one of the tasks for the power penalty approach is proposed for the linear complementarity problem. That is, we will approximate the linear complementarity problem by a nonlinear parabolic PDE in two spatial dimensions with an l k penalty term. In addition, a so-called finite volume method is presented for the numerical solution of the two-dimensional nonlinear partial differential equation. The innovation of this method is in that it combines a finite volume formulation with a fitted approximation. The finite volume method possesses a special feature of the local conservativity of the numerical flux, and is becoming more and more popular. See, for instance, Wang [29] for degenerate parabolic problems, Leveque [30] for hyperbolic problems, and Liu [31] for elliptic problems. Here we aim to present the theoretical valuation of carbon emission derivatives using partial differential equations coupled with numerical methods.
The paper is organized as follows. In Section 2, the pricing PDEs of the CO 2 emission allowance derivatives are obtained, and the final and boundary conditions are prescribed for different contingent claims. Then, a power penalty method is introduced for the American option in Section 3. In Section 4, a so-called fitted volume method is proposed for the discretization of the pricing PDE. Some numerical experiments are performed to illustrate the efficiency and usefulness of the numerical method in Section 5. In addition, some sensitivity analyses are also made to examine the effects of some parameters on the valuation results in Section 6. Finally, some comments are given in Section 7.
The pricing model of CO 2

contingent claims
In order to achieve the pricing equation, which applies to any contingent claims for CO 2 emission allowance spot prices, we assume that the spot price of CO 2 emission allowance S t at time t and the convenience yield δ follow a joint diffusion process specified by (see, for example, [10,32]) dS t ¼ k s mðtÞ Àln ðS t Þ þ 1 k s 1 2 s 2 s þ dm dt À d S t dt þ s s S t dW P s ; ð1Þ where k s and k c denote the speeds of mean-reversion of spot price and convenience yield, σ s and σ c denote the constant volatilities of spot price and convenience yield, respectively, θ c is the long-run mean yield, and μ(t) = ln (the marginal abatement cost) + ξt is the logarithm of the marginal abatement cost C s growing at rate ξ. In addition, dW P s and dW P c are correlated increments to standard Brownian process under the historical measure P, and they have the correlation coefficient ρ sc . From the concept of absence of arbitrage opportunities, if the market prices of spot and convenience yield are given, then we know are Brownian motions under the risk-neutral measure Q, where λ s and λ c are the market prices of spot and convenience yield, respectively. Note that Eq (1) is based on the relationship between the CO 2 emission spot price and the marginal abatement cost given by where X t is governed by an Ornstein-Uhlenbeck process with a zero long run mean, a speed of adjustment k s , and a constant volatility σ s : dX t ¼ Àk s X t dt þ s s dW p s : Eq (3) implies that the natural logarithm of the CO 2 emission allowance spot price deviates around a mean-reversion level, and this level grows with the discount rate ξ and is determined by the marginal cost of abatement C s . Moreover, it is expected to be mean-reversed for convenience yields because of the strong tendency, and short-term random convenience yields converge to their mean values in the long run.
The partial differential equation for the price of carbon derivatives Assume that the price of carbon derivative F(S, δ, t) is a twice continuously differentiable function of S and δ, where we omit the subscript t of S for brevity. Then, we can use Itô's lemma to derive the instantaneous price change as follows: from which the relative change can be obtained via dividing both sides by F For simplicity, we define Following Brennan and Schwartz [33], we construct a portfolio P by investing amounts of x 1 , x 2 , and x 3 in three options of maturities τ 1 , τ 2 , and τ 3 , respectively. Then, the rate of return on this portfolio is given by which is non-stochastic if the portfolio proportions are chosen so that the coefficients of dW P s and dW P c in Eq (6) are zero. That is, To eliminate the arbitrage opportunity, the above return on the portfolio must be risk-free over short time intervals, this is, the return rate is equal to r, the instantaneous risk-free interest rate. As a consequence, the portfolio risk premium is zero: The arbitrage free condition Eq (8) and the two zero risk conditions Eq (7) have a solution only if where λ s and λ c are the market prices of per unit spot and convenience yield, respectively. Subsequently, substituting Eqs (4) and (5) into Eq (9) results in @F @t þ ðk s ðnðtÞ Àln SÞ À dÞS @F @S þ k c ðy c À dÞ @F @d That is, which is the partial differential equation for the contingent claims depending on the evolution of the spot price S and the convenience yield δ. Remark 1. In the case of futures contract, the return rate on the portfolio Eq (6) should be zero [34], such that the Eq (8) changes into x 1 k(τ 1 ) + x 2 k(τ 2 ) + x 3 k(τ 3 ) = 0. Then, the pricing partial differential equation becomes @F @t þ ðk s ðnðtÞ Àln SÞ À d À l s s s ÞS @F @S þ ðk c ðy c À dÞ À l c s c Þ @F @d Remark 2. The same pricing partial differential equation can be also obtained by using twodimensional version of the Feynman-Kac formula.

Boundary and final conditions for different carbon derivatives
It stands to reason that we should prescribe the final conditions for different derivatives to fully define the problem. In addition, to solve the pricing PDE numerically, we also need the boundary conditions. In this paper, we propose three kinds of derivatives, European and American options as well as futures, which are the popular derivatives in any financial market.
Options. We only take the call option's final and boundary conditions into consideration, since the put option is similar. The most common final condition is the vanilla option's payoff given by FðS; d; TÞ ¼ max ð0; S À KÞ; S 2 ð0; S max Þ; ð12Þ where K denotes the exercise price of the option satisfying 0 < K < S max . A second choice is the cash-or-nothing payoff given by where B > 0 is a constant and H denotes the Heaviside function. Obviously, this final condition is a step function, which is zero if S < K and B if S ! K. The price of any derivatives except the futures contract will satisfy Eq (10), and the boundary conditions have to be adjusted according to the specific exercise features of derivatives. In the case of European options, there are four boundaries in the solution domain: S = 0, S = S max , δ = δ min , and δ = δ max . Obviously, the boundary conditions at S = 0 and S = S max are simply taken to be the extension of the final condition: In the case of American options, since they can be exercised at any point in time before maturity, part of the valuation problem involves identifying the optimal exercise policy, or the exercise time that maximizes the option value. The boundary condition at S = S max should be replaced by the following two classic value-matching and smooth-pasting boundary conditions [35]: for the American call option.
To determine the boundary conditions at δ = δ min and δ = δ max , we need to model the pricing Eq (10) again for two particular values δ = δ min and δ = δ max , and solve the two resulting one-dimensional equations. These two approximate boundary conditions will be used to implement our numerical methods.
Futures. In the case of futures, the way to determine the boundary conditions is similar to that for European options. We ignore the discussion about boundary conditions and only propose the final condition in this subsection. As mentioned in [36], when the delivery period is being reached, the futures price is very close or even equal to the spot price, which means that the final condition for futures should be The power penalty approach Since the valuation of American option is a free boundary problem, its exact solution is not available analytically. Therefore, numerical approximation to the solution is normally sought in practice. In fact, we formulate the free boundary problem as a linear complementarity problem, and then develop a power penalty method to solve it. That is, a nonlinear parabolic PDE in two spatial dimensions with an l k penalty term is utilized to approximate the linear complementarity problem. Moreover, in the next section, a fitted finite volume method is proposed to solve the penalized nonlinear equation.

Formulation of the problem into a complementarity problem
Denote the value of the American option with expiry date T by F(S, δ, t) and define where h(S, δ, t) = S(k s (n(t) − ln S) − λ s σ s − δ) and g(δ, t) = k c (θ c − δ) − λ c σ c . Here we define h(0, δ, t) = 0 to keep the continuity of h(S, δ, t) at S = 0. The free boundary S ? (t) divides the domain O = (0, S max ) × (δ min , δ max ) into the continuation region S 1 and the stopping region S 2 . In the continuation region S 1 , F > F(S, δ, T), LF = 0; in the stopping region S 2 , note that S > K, F = S − K, then LF > 0. In a word, the American option value F satisfies the following partial differential complementarity problem: is the payoff function, g 1 , and g 2 are the boundary conditions to be determined in the next section. For the ease of theoretical analysis, we rewrite Eq (18) as the following divergent form: c ¼ r þ k s ðnðtÞ Àln SÞ À d À l s s s À s 2 s À k s À k c : Since the homogeneous Dirichlet boundary conditions are convenient for theoretical discussion, we transform Eqs (19)- (21) to be homogeneous, which is not necessary in computations. To this purpose, we introduce a new function uðS; d; tÞ ¼ e bt ðF 0 À FÞ; ð24Þ is a twice differentiable function satisfying the boundary conditions in Eq (20). Transforming F in Eq (19) into the new function u, we have

Formulation of the complementarity problem into a variational inequality problem
Next we reformulate Eq (25) as a variational inequality problem in an appropriate functional setting. To this end, we first introduce some standard notations. We define O = (0, S max ) × (δ min , δ max ) and Γ as the computational domain and its boundary, respectively. Let Also, we define where u ? (t) is defined by Eq (26). Clearly, K is a convex and closed subset of H 1 0;o ðOÞ. We also define the norm of L p (0, T; H(O)) for any Hilbert space H(O) as follows: k vðÁ; Á; tÞ k L p ð0;T;HðOÞÞ ¼ Hereinafter, v(Á, Á, t) will be written simply as v(t) when there is no confusion. Now, we are in the position to define our variational inequality problem. For this variational inequality problem, we have the following theorem. Theorem 1. Problem 1 is the variational form of the complementarity problem Eq (25).
Theorem 2. Problem 1 has a unique solution.

The power penalty approach
First of all, we introduce the following nonlinear variational inequality to obtain the power penalty method: Seek u l 2 H 1 0;o ðOÞ to satisfy that for all v 2 H 1 0;o ðOÞ, À @u l @t ; v À u l þ Bðu l ; v À u l ; tÞ þ jðvÞ À jðu l Þ ! ðf ; v À u l Þ a:e: in ð0; TÞ; ð33Þ where and [z] + = max{0, z}. According to Lemma 1 and the lower semi-continuity of j, the problem has a unique solution. In addition, j(v) is differentiable, which means that Eq (33) Modelling and Computation in the Valuation of Carbon Derivatives Note that the penalized variational equation of variational inequality Eq (27) is just Eqs (35)- (36), and the penalized equation, which approximates Eq (25), can be written as: with the following boundary and final conditions: It has been proved in several papers [37][38][39][40] that the solution of Eq (37) converges to that of Eq (19) at the rate O(λ − k/2 ) in a Sobolev norm as λ ! + 1.

The fitted finite volume method
As the valuation models for futures and European options are the special forms of Eq (37), the so-called fitted finite volume method is presented only for the American option model here, and it can be also applied to the other two models. The innovation of this method is that it combines two existing techniques, a finite volume method and a fitted approximation, together, in which the flux of a given function is approximated by a constant locally, yielding a locally nonlinear approximation to the function. In what follows, we will develop the method for our two-dimensional nonlinear partial differential equation with a penalty term.

Boundary conditions
First of all, how to determine the boundary condition functions g 1 (S, t) and g 2 (S, t) is discussed below.
1. The boundary condition g 1 (S, t) on the boundary δ = δ min is determined by solving the following parabolic partial differential equation: FðS; TÞ ¼ max ðS À K; 0Þ: 2. The boundary condition g 2 (S, t) on the boundary δ = δ max is determined by solving the following initial-boundary problem: FðS; TÞ ¼ max ðS À K; 0Þ: The fitted finite volume method From Eq (24) we can easily find that the problem Eqs (37)- (38) can be rewritten as with the boundary and the terminal conditions Eqs (20)- (21). In what follows, we let I S = (0, S max ) and I δ = (δ min , δ max ), and divide the intervals I S and I δ into N S and N δ sub-intervals, respectively: This defines a mesh on I S × I δ with all mesh lines perpendicular to one of the axes. Next we define another partition of I S × I δ by letting for i = 1, 2, Á Á Á, N S − 1 and j = 1, 2, Á Á Á, N δ − 1, and S À 1 Applying the mid-point quadrature rule to the above equation, we have for i = 1, 2, Á Á Á, N S − 1, j = 1, 2, Á Á Á, N δ − 1, where R i;j ¼ ðS iþ 1 2 À S iÀ 1 2 Þ Â ðd jþ 1 2 À d jÀ 1 2 Þ, c i;j ¼ cðS i ; d j ; tÞ, F i,j = F(S i , δ j , t), and F Ã i;j ¼ F Ã ðS i ; d j ; tÞ: Now we consider how to approximate the second term in Eq (42). First of all, it follows from the definition of flux ArF þ b F and integration by parts that Z where l denote the unit vector outward-normal to @R i,j . Next we deal with Eq (43) term by term. In fact, for the first term we can approximate the integral by a constant, i.e, Z S iþ 1 2 Clearly, we now need to derive approximations of the ðArF þ bFÞ Á l defined above at the mid-point, ðS iþ 1 2 ; d j Þ, on the interval I S i for any i = 0, 1, Á Á Á, N S − 1. Next we discuss it in detail. Case 1: For i ! 1. According to Eq (23) we have Following the idea in [29], we approximate the term a 11 F S + b 1 F by solving the following twopoint boundary value problem: where a ¼ 1 2 s 2 s and b ¼ k s ðnðtÞ À lnSÞ À d À l s s s À s 2 s , b iþ 1 2 ;j ¼ bðS iþ 1 2 ; d j Þ, and F i,j and F i + 1,j denote the values of F at (S i , δ j ) and (S i + 1 , δ j ), respectively. Evidently, from Eq (45a) we can obtain where C 1 is an arbitrary constant and can be determined by the boundary condition Eq (45b) as follows [29,41]: a . Then, where d ¼ 1 2 r sc s s s c . Additionally, the derivative F δ can be approximated by a forward difference Then, we have where d i,j = d(S i , δ j ). In the same way, we can also approximate the second term in Eq (43) as follows: Case 2: For i = 0. As the Eq (45a) is degenerate on (0, S 1 ), we should rethink the problem Eqs (45a)-(45b) on (0, S 1 ) below: where b1 2 ;j ¼ bðS1 2 ; d j Þ and C 2 is an unknown constant to be determined next. After integrating both sides of above equation, we can obtain Thus, we have In the approximations of last two terms of Eq (43), we do not need to consider the case 2 as before, since δ 0 6 ¼ 0. We apply the similar method to the above and get the following results: for j = 0, 1, Á Á Á, N δ − 1, where a i;j ¼ b i;jþ 1 2 a j , a ¼ 1 2 s 2 c , b ¼ k c ðy c À dÞ À l c s c À 1 2 r sc s s s c , and d i;j ¼ 1 2 r sc s s s c S i . Hence, we obtain the following equations by combining Eqs (43), (47), and (48) with Eqs (50)-(52): À @F i;j @t R i;j þ e i;j iÀ1;j F iÀ1;j þ e i;j i;jÀ1 F i;jÀ1 þ e i;j i;j F i;j þ e i;j i;jþ1 F i;jþ1 for i = 1, 2, Á Á Á, N S − 1 and j = 1, 2, Á Á Á, N δ − 1, where for j = 1, 2, Á Á Á, N δ − 1, and ; e i;j i;jÀ1 ¼ À b i;jÀ 1 2 e a i;jÀ1 d jÀ1 h S i e a i;jÀ1 d j À e a i;jÀ1 d jÀ1 ; ð57Þ for i = 2, 3, Á Á Á, N S − 1, j = 1, 2, Á Á Á, N δ − 1, and e i;j m;n ¼ 0 if m 6 ¼ i − 1, i, i + 1 and n 6 ¼ j − 1, j, j + 1. It can be easily seen that Eq (53) is an (N S − 1) 2 × (N δ − 1) 2 linear system of equations for F ¼ ðF 1;1 ; Á Á Á ; F 1;N d À1 ; F 2;1 ; Á Á Á ; F 2;N d À1 ; Á Á Á ; F N S À1;1 ; F N S À1;2 ; Á Á Á ; F N S À1;N d À1 Þ > : Note that for i = 1, 2, Á Á Á, N S and j = 1, 2, Á Á Á, N δ , F 0, j (t), F i,0 (t), F 0, N δ (t) and F N S ,0 (t) are equal to the given boundary conditions. Obviously, the coefficient matrix of Eq (53) is penta-diagonal. Let E i;j ¼ ð0; Á Á Á ; 0; e i;j iÀ1;j ; 0; Á Á Á ; 0; e i;j i;jÀ1 ; e i;j i;j ; e i;j i;jþ1 ; 0; Á Á Á ; 0; e i;j iþ1;j ; 0; Á Á Á ; 0Þ for i = 1, 2, Á Á Á, N S − 1 and j = 1, 2, Á Á Á, N δ − 1. Now we are in the position to discuss the time-discretization of the system Eq (53). To this purpose, we first rewrite Eq (53) as Then, we select M − 1 points numbered from t 1 to t M − 1 between 0 and T and let T = t 0 , t M = 0 to form a partition: T = t 0 > t 1 > Á Á Á > t M = 0. Thus, the full discrete form of Eq (60) can be obtained by applying the two-level implicit time-stepping method with a splitting parameter y 2 ½ 1 2 ; 1 to it as follows: where for

The solution of the discrete system
The standard Newton method is employed to solve the nonlinear discrete system Eq (62). Eq (61) clearly indicates that p 0 ðV m i;j Þ ! 1 as F Ã i;j À F i;j ! 0 þ when k > 1. So, we need to smooth out pðF m i;j Þ by redefining it as follows: where k > 0, n is a positive integer, and 0 < ( 1 is a transition parameter. Applying Newton method to Eq (62), we can get for q = 1, 2, Á Á Á, where ω 0 is a given initial guess, J D (ω) denotes the Jacobian matrix of the column vector D(ω), and γ 2 (0, 1] is a damping parameter used to accelerate convergence. Then, we choose as the solution of Eq (65). Remark 4. For the European options and futures, the system Eq (62) degenerates to a linear system and can be solved using normal methods.
What is more, we have the following theorem. Proof. From the definition of D(F) in Eq (63), it is easy to see that its Jacobian is the following diagonal matrix: From Eq (64) we know that p 0 ðF m i;j Þ ! 0 for all i = 1, Á Á Á, N S − 1 and j = 1, Á Á Á, N δ − 1. Thus, to show that the system matrix of Eq (65) is an M-matrix, it suffices to show that θE m + 1 + G m is an M-matrix.
First, we note that e i;j m;n 0 for all m 6 ¼ i, for any i and j, and for any α = b/a and any a ¼ b= a. This is because the function S α is increasing when b > 0 and decreasing when b < 0, and the function e ad is increasing when b > 0 and decreasing when b < 0. Eq (66) also holds when b iþ 1 2 ;j ! 0, b i;jþ 1 2 ! 0. Furthermore, from Eqs (54)-(56) we know that when c i;j ! 0, for all i = 1, Á Á Á, N S − 1, j = 1, Á Á Á, N δ − 1, there holds jðe i;j p;q Þ mþ1 j þ c mþ1 i;j R i;j : Therefore, E m + 1 is a diagonally dominant with respect to its columns. Hence, from the above analysis, we see that for all admissible i, j, E m + 1 is a diagonally dominant matrix with positive diagonal elements and non-positive off-diagonal elements. This implies that E m + 1 is an Mmatrix. Second, G m of the system matrix Eq (65) is a diagonal matrix with positive diagonal entries. In fact, when jΔt m j is sufficiently small, we have which demonstrates that θE m + 1 + G m is an M-matrix.
Here we emphasize that Theorem 3 implies that the fully discrete system Eq (65) satisfies the discrete maximum principle and the discretization is monotone, such that Eq (65) has a unique solution.

Numerical results
In this section, some numerical results are presented to demonstrate the efficiency and the usefulness of the numerical method proposed in the above. Furthermore, the varieties of derivatives prices with respect to the spot price, convenience and time, are examined. These are of interest to market participants. To solve the pricing problem numerically, we divide (0, S max ), (δ min , δ max ), and (0, T) uniformly into 50, 50, and 50 sub-intervals, respectively. The boundary conditions are given by Eqs (14), (39), and (40), and the numerical values of these boundary conditions determined by 1D initial-boundary problems are plotted in Fig 1. By means of these initial and boundary conditions, the European option values for Test 1 are computed and the intersection surfaces of different time points are plotted in Fig 2. We show some option values at some special points in Table 1.
Next the convergence rates of the discretization method is gauged. To this end, we define three discrete norms k Á k 1,h S , k Á k 1,h δ , and k Á k 0,h as follows [42]: from which we can define the following weighted discrete H 1 -norm:  In addition, the ratio is defined as: In Test 1, we employ the numerical solution on the mesh with N S = 128 = N δ and M = 128 as the "exact" solution F, and list the errors in the discrete L 1 -norm and the weighted discrete H 1 -norm at the final time step t = 0 for four consecutive meshes in Table 2. Moreover, the linear regression is used to show that these data obey the basic error estimates as follows: k F À F h k 1 % 0:3136h 1:0597 and k F À F h k H 1 % 2:0922h 1:6493 : Note that these results are reasonable because of the non-smoothness of the solution due to the pay-off function Eq (13). From Figs 1, 2, and Table 1 we can conclude that 1. The European cash-or-nothing option values are higher when the option comes to maturity. This is acceptable in that the European option can be only exercised at maturity, and one can not get more benefits when holding a European option that is not closed to maturity, which is different from the American option.
2. The greater convenience yield δ, the lower the option value. This can be explained as follows: since the convenience yield δ measures the benefits from holding the carbon emission permits rather than the derivatives contract, people will hold more carbon emissions spots when facing a higher convenience yield, which leads to the demand reduction for derivatives and the option price down.
3. As the value of carbon emission spots rises, the value of call option also rises. This is because the holder of call option will get more chances for benefits from the derivatives. This result is similar to the general call options. To solve the American call option with the above parameters in Test 2, we divide (0, S max ), (δ min , δ max ), and (0, T) uniformly into 30, 30, and 30 sub-intervals, respectively. The final and boundary conditions are given by Eqs (20), (21), (39), and (40), and the numerical values of these boundary conditions determined by the 1D initial-value problems are plotted in Fig 3. According to these initial and boundary conditions, the American call option problem for Test 2 is solved and the intersection surfaces of different time points are depicted in Fig 4. We also show some option values at some special points in Table 3. From Figs 3, 4 and Table 3, we have the conclusions as follows: 1. The American vanilla call option values are lower when the option comes to maturity. An American option can be exercised at any time between the date of purchase and the expiration date. As the option comes to maturity, the benefits from the early exercise of option are reduced, and then the option prices decrease. This is the same as the popular financial option contracts.
2. Unlike European options, for American options the greater convenience yield δ does not mean a lower option value (see the last two lines of Table 3). As we all have known, the value of American option contains an early exercise premium. Therefore, although the greater convenience yield reduces the demand for derivatives and decreases the option price, the early exercise premium may also increase the option price oppositely. So, there is no certain relationship between the convenience yield and the option value.
3. As the value of carbon emission spots rises, the value of the call option also goes up. This is because the holder of the call option can get more chances to benefit from it. This result is similar to general call options.
Financial engineers could pay more attention to the optimal exercise boundary compared with the option value in the case of American option. We determine the boundary S ? (t) by Eq (16), and the results are plotted in Fig 5. The blue line is the optimal exercise boundary, the left of the blue line is the continuation region, and the right of the blue line is the stopping region. We can clearly see that as the convenience yield increases, the stopping region becomes smaller. Note that our pricing model is different from the classical Black-Scholes model in details, and this result should be reasonable.
Test 3: The futures with the final condition Eq (17).    To solve the futures with the above parameters in Test 3, we divide (0, S max ), (δ min , δ max ), and (0, T) uniformly into 50, 50, and 50 sub-intervals, respectively. The boundary conditions are given by Eqs (14), (39), and (40), and the numerical values of these boundary conditions determined by the 1D initial-value problems are plotted in Fig 6. According to these initial and boundary conditions, the futures problem for Test 3 is solved, and the intersection surfaces of different time points are depicted in Fig 7. Some futures values at some special points are shown in Table 4. From Figs 6, 7 and Table 4, we have the conclusions as follows: 1. The futures prices are lower when the futures contract comes to the maturity. This is because the futures prices converge to the spot prices when the futures contract comes to the maturity. With time approaching to the maturity date, the uncertainty, which should increase the futures prices, is reducing.
2. The greater the convenience yield δ, the lower the futures value. This can be explained as follows: since the convenience yield δ measures the benefits from holding the carbon emission permits rather than the derivatives contract, people will hold more carbon emission spots when facing a higher convenience yield, which leads to the demand reduction for derivatives and the futures price down.
3. As the value of carbon emission spots goes up, the value of futures also rises. This is because the holder of futures will get more chances to benefit from the derivatives. This result is similar to general futures.

Discussions
In this section, we examine the sensitivity of above results to the parameters variations. There are four parameters which are likely to have an impact on the price of emission allowance derivatives. Discussion 1. We examine the effects of some parameters on the European cash-or-nothing option prices in Fig 8, where we fix t = 0 and δ = − 1.
From Fig 8 we can see that a. The option price increases with the higher growing rate of marginal abatement costs. The higher growing rate of marginal abatement costs will make the abatement costs change quickly and the enterprises will live in a more uncertain world when they must choose investing on clean technology or trading the emission permits. To eliminate this uncertainty, enterprises will buy more derivatives such as options, which leads to the option value increasing. b. The effect of the marginal abatement costs is negligible. The purpose to trade the emission permits is just to hedge the risk of marginal abatement costs, and if marginal abatement costs increase, people will hold more emission permit spots, and the spots price will increase and vice versa. So, it is natural that the option value does not change following the movement of marginal abatement costs.
c. Increasing the speed of mean-reversion k s reduces the option prices. This is because emission permit prices have a stronger tendency towards the mean value, the price risk of emission permits decreases, and then people need not to hold many derivatives to hedge the risk. The decrease of demand will reduce the option prices.  d. The effect of the volatility σ s is irregular compared with the classical Black-Scholes model. In the Black-Scholes world, the higher volatility of underlying assets increases the risk of investment, and the derivatives price will increase. In our model, since the stochastic convenience yield is considered, the effect of volatility on option prices may be counteracted by the stochastic convenience yield. So, there is no certain relationship between the volatility and the option price when there is a stochastic convenience yield.
Note that the above sensitivity analysis results are also satisfied by the following American option and futures, as they are all the derivatives used to hedge the price risk.
Additionally, the effects of stochasticity of convenience yield on the valuations of carbon European options are also examined in Fig 9, where t = 0 is fixed. Note that the solid line states the deterministic convenience yield simulation results and the dash line represents the stochastic ones.  c. Increasing the speed of mean-reversion k s reduces the option prices.
d. The effect of the volatility σ s is irregular compared with the classical Black-Scholes model.
Similarly, the effects of stochasticity of convenience yield on the valuation of carbon American option are also examined in Fig 11, where t = 0 is fixed. Note that the solid line states the deterministic convenience yield simulation results and the dash line represents the stochastic ones.  c. Increasing the speed of mean-reversion k s reduces the futures price, because the emission permits price trends to the mean value more quickly, and the uncertainty of emission permits price would become smaller.
d. There is no certain effect of the volatility on futures prices. Once again, the effects of stochasticity of convenience yield on the valuations of carbon futures are examined in Fig 13 with the fixed t = 0. Note that the solid line corresponds to the deterministic convenience yield simulation results and the dash line represents the stochastic ones.
From the above discussions, we can summarize some experiences about risk management. Firstly, as the Kyoto Protocol came into force in 2005, several countries are expanding their investment scales on the new technologies to reduce carbon emissions, and thus the emission allowance derivatives are also becoming a risk hedging tool for market participants. Secondly, a firm may not be scarce to the emission allowance derivatives when it can adjust its marginal abatement costs well. Finally, the impact of volatility of carbon emission allowances price is difficult to measure when a stochastic convenience yield is involved.

Concluding remarks
This paper presents a methodology for modelling and computing the valuation of carbon derivatives with stochastic convenience yields. The principle of absence of arbitrage opportunities  and the stochastic calculus are used to develop the mathematical model, the partial differential equation, satisfied by carbon derivatives. For American options, we formulate the pricing problem to a linear parabolic variational inequality, and develop a power penalty method to solve it. Then, a so-called fitted finite volume method is designed to solve the nonlinear partial differential equation resulting from the power penalty method, which governs the futures, European and American option valuation.
Also, we discuss the effects of the stochastic convenience yield on the prices of carbon derivatives in details. First of all, a greater convenience yield means a lower European option price and futures price while it does not mean a lower American option price. Therefore, the property of optimal exercise boundary makes the American option quite different from the European option. Furthermore, from the sensitivity analysis we can see that (a) the derivatives price increases with the higher growing rate of marginal abatement costs; (b) the effect of the marginal abatement costs is negligible; (c) increasing the speed of mean-reversion k s reduces the The carbon market participants, such as investors, hedgers and arbitragers, may get help from our theoretical results, and they can work out better hedging strategies, adjust portfolio structures and strengthen their capabilities to manage risks. We expect that our work can emphasize the importance of the convenience yield concept for the emission allowance pricing. At the same time, we also anticipate that our derivatives pricing methodology from the perspective of partial differential equations combined with numerical methods can make a few contributions to the development of carbon market.

Author Contributions
Conceived and designed the experiments: SZ XW. Performed the experiments: SZ XW. Analyzed the data: SZ XW. Contributed reagents/materials/analysis tools: SZ XW. Wrote the paper: SZ XW.