MINIMIZATION OF CARBON ABATEMENT COST: MODELING, ANALYSIS AND SIMULATION

. In this paper, we consider a problem of minimizing the carbon abatement cost of a country. Two models are built within the stochastic op- timal control framework based on two types of abatement policies. The corresponding HJB equations are deduced, and the existence and uniqueness of their classical solutions are established by PDE methods. Using parameters in the models obtained from real data, we carried out numerical simulations via semi-implicit method. Then we discussed the properties of the optimal policies and minimal costs. Our results suggest that a country needs to keep a relatively low economy and population growth rate and keep a stable economy in order to reduce the total carbon abatement cost. In the long run, it’s better for a country to seek for more eﬃcient carbon abatement techniques and an environmentally friendly way of economic development.


1.
Introduction. Global warming, as a result of greenhouse gas (GHGs) emissions, is causing severe consequences all around the world. Researchers have shown that the long term cost due to global warming will be significantly higher than that of controlling it by reducing GHG emissions(McCarthy et al. [35]). In 1997, international communities adopted the Kyoto Protocol to encourage countries participate in GHG emission reduction, or carbon reduction. The Kyoto Protocol has set emission limits on Annex B nations, and offered three mechanisms, including international emission trading (IET), for them to comply with their commitments (Carmona et al [6]; Fehr and Hinz [19]). In Jun. 2000, the European Commission (EC) launched European Climate Change Program (ECCP) to implement the Kyoto Protocol in the European Union (EU)(Belauar et al [4]; Cetin and Verschuere [8]). Then ECCP implemented the European Union Emission Trading Scheme (EU ETS), which began operating since Jan. 1, 2005(Ellerman and Buchner [15]), and is now the largest carbon market in the world(Hepburn [23]).
The EU ETS runs a cap-and-trade scheme by imposing an upper bound on the total emission of every country and allowing these countries to trade their emission allowances (Belauar et al [4]; Carmona et al [7]; European Commission [12]). Countries included in the EU ETS need to surrender allowances equivalent to their verified CO 2 emissions every calendar year, or they have to pay penalties for the excess emissions(European Commission [11,12]). This leads to a bilateral transaction called "carbon trading".
With these subjects on carbon trading becoming popular, researchers have carried out studies on them from many aspects. Some of these works carried out by Cetin and Verschuere([8]), Carmona, Fehr and Hinz( [6]), Carmona, Fehr, Hinz and Porchet( [7]), etc, directly analyzed the price of traded emission permits and the behaviour of the market. Meanwhile, some works, like Zagheni and Billari( [47]), Belauar, Fahim and Touzi( [4]), etc, considered the choice of a country or company to reduce its carbon emission at a reasonable cost in the carbon market. Zagheni and Billari( [47]) used the Stochastic Impacts by Regression on Population, Affluence and Technology (STIRPAT) equation proposed by Dietz and Rosa([14]) and built a cost valuation model to assess the economic costs of a country related to carbon emission reduction by option pricing method. Belauar, Fahim, Touzi( [4]) considered an installation with the obligation to reduce its carbon emission. The installation aims at obtaining the maximal wealth from their products and carbon emission trading by choosing an optimal production process. The carbon emission process of the factory is modeled as a deterministic function of the production process, while the carbon emission market and the value of carbon emission contracts are modeled by stochastic processes.
In this paper, we consider the problem of carbon emission abatement from the aspect of contingent claims. We use a stochastic process to describe the carbon emission of a country, and find out the optimal carbon abatement policy to minimize the total cost in reducing carbon emission. In this paper, the problem is simplified to a Merton optimal consumption problem and is solved by Hamilton-Jacobi-Bellman (HJB) equations. As far as we know, our earlier work(Yang and Liang [46]) is the first model to solve the problem in this way. Compared with the preliminary model in [46], the model in this paper is more general and more realistic.
HJB equations are nonlinear partial differential equations (PDEs) arising from dynamic programming principle in optimization problems. It is common to seek continuous viscosity solutions (see Fleming and Soner [21] for more details). In our model, the HJB equations happen to be semilinear PDEs with bounded coefficients, which makes it possible to seek for classical solutions. In contrast to [21], our equations are on unbounded domains, the solutions themselves are unbounded. The initial condition contains a Lipshitz cusp. In Model 2, the first order derivative of the coefficient also becomes unbounded. We shall extend classical PDE results to these equations. In this paper, the existence and uniqueness of the classical solutions are proved to show that the model is rational, and properties of them are shown by numerical results and discussed.
The rest of the paper is organized as follows. In section 2, the model of minimizing the total cost of carbon emission reduction and the penalties due to extra carbon emission is built, and the HJB equation is deduced. In section 3, we study the problem of reducing carbon emission growth rate and prove the existence, uniqueness and regularity of the solution. In section 4, we study the problem of reducing carbon emission amount and establish the existence and uniqueness of the solution. In section 5, we describe the data used to obtain the parameters in the models. In section 6, we present some simulation results. And finally in section 7, we summarize what our model can predict.

2.
Modeling. Let (Ω, F , {F t } t≥0 , P) be a filtered probability space. In this paper we consider a country which is obliged to reduce its domestic carbon emission. From the IPAT model(Dietz and Rosa [14]; Zagheni and Billari [47]), The contributors of the carbon emission (I) of an area are population (P ), Affluence (Economic good/Population, A) and technology (Pollutant/Economic good, T ). The model is often used to solve T with given data on I, P and A(Dietz and Rosa [14]). Then Dietz and Rosa([14]) reformulated the model into a nonlinear one and made it possible to evaluate the weights of different contributors by statistical methods. The contributors to carbon emission in this model are P and A, and T is modeled as a residual term. Then Zagheni and Billari( [47]) replaced the factor A (represented by per capita GDP) by Y = AP , i. e. GDP, to build a stochastic differential equation (SDE for short) to describe the movement of carbon emission. This modified STIRPAT model was later used in [46], and will be the model used in this paper as well.
In the modified STIRPAT model, the carbon emission of a country is determined by the economy and the population, and GDP can be used to represent the economy of the country. Denote the GDP as Y t and assume Y t satisfies geometric Brownian motion: where µ and σ are both positive constants, and W t is an F t -adapted standard Brownian motion.

Differentiated ln(GDP)
Year Differentiated ln(GDP)   Figure 1 (a) shows a clear trend in the time series; clearly the time series is not stationary. This is confirmed by the ACF figure (lower part of 1 (a)). In this figure, the ACF of the time series decreases slowly to zero instead of decreasing exponentially, which is a property of stationary time series. However, the log-return series of the annual GDP (differentiated ln(GDP)) is more stationary. This can be confirmed by the plot of the series and ACFs in Figure 1 (b). Thus, a Brownian motion is a good approximation of the log-return. Then the annual GDP follows the geometric Brownian motion. Similar practices of using geometric Brownian motion to describe GDP movement can be found in Schinckus [40]; Kruse et al [33]; Chamon and Mauro [9]; Zagheni and Billari [47].
Assume that the population P t of the country satisfies the Logistic model (see [43] for detailed descriptions). Assume that the population of an area at time t is P t , the population at time t = 0 is P 0 , the carrying capacity of the population of the area is P m , and the intrinsic population growth rate is ρ. In a Logistic model, the actual population growth rate should be By solving the above ordinary differential equation (ODE) we can obtain the expression of P t as follows: It can be found that P t satisfies the following equation: In a modified STIRPAT model in [46], the relationship between the carbon emission I t of the country, the GDP and the population is as follows: where a 1 , a 2 are assumed to be constants representing the contributions of GDP and population in affecting the carbon emission respectively. Using (2.1)(2.2) in (2.3), we find that the evolution of the carbon emission satisfies the following SDE: Two kinds of policies are applied to reduce carbon emission. The first kind is to reduce the growth rate of emission, while the second one is to reduce the amount of emission directly. In this paper, two models will be established based on these two types. For decision makers to observe the annual carbon emission growth rate and assess the current carbon abatement policies from a macroscopic aspect, the first model is better. However, the second model is more practical by making the modeling of carbon abatement costs much easier. There are more mathematical difficulties in this model.
Model 1: Reduce emission growth rate Assume the reduced increasing rate of carbon emission is C t , and equation (2.4) is transformed into: In order to keep the cost of reducing carbon emission from becoming too high, an upper bound is placed on C t . TakeC as the upper bound, which is a positive constant, then the (admissible) control policy C t is assumed to be an {F t } t≥0adapted process satisfying 0 ≤ C t ≤C for all t ≥ 0. And the set of all admissible controls is denoted by C.
The total cost of reducing carbon emission comes from two aspects: (1) the cost from lowering carbon emission growth in the whole period 0 ≤ t ≤ T when the control is in place, and (2) the penalty at termination T if the carbon emission quota is exceeded. The first kind of cost is calculated continuously with a cost rate assumed to be g 1 (C t ), meaning that from time t to time t + dt the cost is g 1 (C t )dt. It is reasonable to assume that g 1 (C) is a non-negative, increasing function of C. The implications of this assumption is simple: first, the cost to reduce carbon emission growth rate is non-negative; second, the more emission to reduce, the higher cost would be needed.
The risk-free interest rate is assumed to be constant r > 0. Then the present value of the total cost of reducing emission from time t(< T ) to termination T is T t e −r(s−t) g 1 (C s )ds.
At termination T , if the carbon emission of the country exceeds a constant thresh-oldĪ > 0 set at time t = 0, the country has to pay a 0 (positive constant) for every ton of extra emission. So the cost of carbon emission at time T should be a 0 (I T −Ī) + . And the conditional expectation at time t of the country's total cost of reducing carbon emission should be The goal of a country is to choose the optimal policy of reducing carbon emission growth to make sure that the total cost of reducing carbon emission together with the penalties at termination attains its minimum, i.e. to choose the optimal C * t that minimizes J 1 (I, t; C). Here we define as the minimal cost. If Φ(I, t) is continuous, then it is a viscosity solution to the following HJB equation (c. f. [21]): (2.6) Model 2: Reduce emission amount Assume the reduced amount of carbon emission at time t is A t . Then the controlled carbon emission process (2.4) becomes (2.7) Assume the upper bound of the policy is a positive constantĀ. Define admissible policies as above, and denote the set of admissible policies by A. Assume the cost rate of reducing emission is g 2 (A), which is a non-negative, increasing function of A, and the total cost functional is As in the first model, we seek to minimize the total cost: Similarly, Ψ(I, t) is a viscosity solution of the following equation 3. Model 1: reduce emission growth rate. Applying logarithmic transformation to the HJB equation derived from Model 1 by taking x = ln I, τ = T − t, still denoting τ by t, we derive the following from equation (2.6): and Then the above equation can be rewritten as 3.1. Existence. Recall that in (2.1)-(2.5), µ, σ, r, a 0 are positive constants, a 1 , a 2 are constants, f (t) is defined as where ρ, P m , P 0 are positive constants, C t is the control policy appeared in the formula (2.5). Recall also we have assumed that g 1 (x) is a non-negative, continuous, increasing function of x. Here is our existence result: There exists a solution Φ ∈ C 2+α,1+α/2 (Q ∞,T ) ∩ C (Q ∞,T ) to equation (3.3). Furthermore, the solution Φ belongs to M 1 , where Proof. We will carry out our proof using the following procedures: 1), Build a sequence of solutions on bounded domains with approximate smooth boundary and initial conditions; 2), Obtain uniform estimates of these solutions by a priori estimate; 3), Take a limit to obtain the solution on an unbounded region with the non-smooth initial condition. The proof is as follows: Step 1. Solution on bounded domain Define Ω R = (−R, R), Q RT = Ω R ×(0, T ]. Without loss of generality, assume that R > lnĪ. Consider the following problem: We approximate the initial and boundary conditions so that the second-order compatibility conditions are satisfied: are defined as follows: These functions satisfy where is a small positive constant. Problem (3.5) admits a classical solution (see [10,Thm 6.2]). The solution to (3.5 In particular, h(t) and h (t) are bounded. Define Since can be taken small enough so that < a 0Ī , it is easy to verify that W U and W D satisfy the following inequalities in Q RT : On the parabolic boundary ∂ p Q RT , W U and W D satisfy It can be proved that, for any p, q ∈ R, (3.8), by the strong maximum principle ([25, Thm 3.7]), the non-negative maximum of Φ R − W U can only be attained on ∂ p Q RT . In particular, Thus we established the following inequality: We can rewrite (3.3) as: we apply [10, Thm 9.6.2] to conclude that there exist at least one solution Φ R ∈ C 2+α,1+α/2 (Q RT ) to the initial-boundary problem (3.5).
Step 2. Uniform C 2+α,1+α/2 -estimate For any given L > 0, consider Q LT and the sequence [25,Thm 3.10], we find that the following estimate holds on Q where C depends only on p, Q K , Q LT . In particular, by the embedding theorem([25, Thm 3.14]), for p > 3 and α = 1 − 3/p, where C depends only on p, Ω K .
Using Φ Rx ∈ C α,α/2 (Q K ), and the Lipschitz continuity of G 1 (p), we derive that . In order to obtain higher order derivative estimates, we shall employ Schauder's interior estimates. This means we have to shrink our domain from Q k to Q K/2 . Notice that this will not be a problem since K is arbitrary.

Step 3. Solution on unbounded domain
This is the typical diagonalization process. Take K = 1.
Repeat this process, and we obtain subsequences We choose the diagonal subsequence {Φ R kk } k≥1 and denote the limit of this sequence by Φ 1 , then Φ 1 satisfies equation (3.4) We repeat the process to obtain subsequences We choose again the diagonal subsequence Step 4. Φ satisfies the initial condition We shall use the barrier method. Let We can verify that U (x, t) is the classic solution to the following equation (see Jiang [28] page 74-89 for detailed information): i. e. U is the expected total cost on extra carbon emission at terminal time if no carbon abatement measures are taken. From above it is clear that U ∈ C 2+α,1+α/2 (Q ∞,T ), then U is Lipschitz continuous in Q ∞,T . U and ∂U ∂x are both positive and bounded in Q ∞,T : A direct computation shows Similarly, R is taken large enough such that where , e (a1µ+a2ρ−r)T , e aT + 1 , For any L > 0, we consider Φ R with R > L and U defined above, and we use 1 Ω to denote the characteristic function on Ω, and the parameters are chosen as follows: Notice that W

and its minimum is given by
Then we have in Q RT :

XIAOLI YANG, JIN LIANG AND BEI HU
Furthermore, when R is large enough such that 1 2 e n1(R−L) ≥ n 1 (R − L) + 1, the following inequalities hold on ∂ p Q RT : From the above inequalities, and the fact that the x-derivative of W 1 1 is continuous, it can be deduced that Φ R −U −W 1 1 is a weak sub-solution (A weak sub-solution, when compared with a weak solution or a solution, satisfies a comparison principle. see [25,Def 3.2] for more information) to the following equation: By the weak maximum principle ([25, Similarly, it can be verified that Assuming that M 4 ≤ e R , then it can be proved that Let → 0, R → +∞, and it is obtained that Φ(x, 0) = ψ 0 (f 0 (x)) for x ∈ Ω L . Let L → +∞ and we have Φ(x, 0) = ψ 0 (f 0 (x)) for all x ∈ R.
The above proof leads to the following corollary: Corollary 1. If g 1 (0) = 0, then the solution Φ to equation (3.3) satisfies the following inequality: Recall that U represents the total cost of carbon emission without control. And it is natural to assume that g 1 (0) = 0, i. e., if no measures are taken to reduce carbon emission, the carbon abatement cost should be zero. Then the corollary indicates that the total minimum cost of reducing carbon emission is no greater than that without any carbon abatement policies.

Uniqueness. Define
We have shown that a classical solution exists in the class M 1 . Now we show that the solution is actually unique in the same class. Proof. Assume that there exist two solutions Φ 1 , Φ 2 ∈ M 1 to equation (3.3), then there exist positive constants C 1 , In addition, Φ 1 − Φ 2 satisfies the following equation: In Q RT , W 1 2 satisfies We apply strong maximum principle ( [25,Thm 3.7] Let R → +∞ and we have (Φ 1 − Φ 2 )(x, t) = 0. Therefore Φ 1 − Φ 2 ≡ 0 on Q ∞,T , i.e. the solution to the problem in M 1 is unique. 4. Model 2: reduce emission amount. We apply the logarithmic transformation to equation (2.9) by taking x = ln I, τ = T − t (still denoted by t) to obtain the following equation: where where ρ, P m , P 0 are positive constants, A t is the control policy appeared in the formula (2.8). g 2 (x) is a non-negative, increasing function of x. We further assume that g 2 (0) = 0. Below is our existence result.
Proof. With an unbounded coefficient of the first-order derivative in the nonlinear term, we cannot apply an existence theorem directly to build a solution on a bounded domain with smooth definite conditions, as in the proof of Theorem 3.1. We shall use the contraction mapping principle. Thus the main steps of establishing the existence result are: 1), Consider the solution to an auxiliary equation, and use contraction mapping principle to obtain a solution on a bounded region. 2), Obtain uniform interior-initial estimates of a sequence of solutions on bounded domain, and establish the solution to equation (4.1).

Step 1. Auxiliary solution F R on bounded domain
As before, define Ω R (−R, R), Q RT Ω R × (0, T ]. Without loss of generality, assume that R > lnĪ. We consider the following problem and denote the solution to this problem by Ψ R : Take W U = a 0 e at+x , W D = 0, and a as defined in (3.7). Similar to that in section 3.1, we can prove that W D ≤ Ψ R ≤ W U in Q RT by applying the strong maximum principle ( [25,Thm 3.7]). This leads to the following inequality: (3.10). In Q RT , F R is the solution to the following problem: with 0 < τ < T /2 to be determined. For any v ∈ B, we define F v as solution to the following equation: Since G 2 (x, 0) = 0, we have Similar to (3.8), it can be proved that By applying [25, Thm 3.11] and the embedding theorem ([25, Thm 3.14]), we have the following estimate for p > 3 and α = 1 − 3/p: then it suffices to prove that T is a contraction. First of all, we have the following proposition: Proof. Let F v1 = T v 1 , F v2 = T v 2 , then F v1 , F v2 are solutions to the following equations: x ∈ Ω R , i = 1, 2.
Step 2. Regularity of F R For any L > 0, we choose R > D = 2L and denote by F R the solution to equation (4.4). From (4.3), it can be shown that , (x, t) ∈ Q DT . For any α ∈ (0, 1), let β = 1+α 2 and choose ξ(x) ∈ C 2+α (Ω D ) such that ξ satisfies the following conditions: By applying [25,Thm 3.11], we derive the following estimate: where C depends only on p, Q DT . With the boundedness of H R , ξ, ξ , ξ , G 2 (x, U x ), there exists a constant C > 0 depending only on p, Q DT such that
Similar to Corollary 1, it can be proved by barrier method that where U is defined in (3.10). This means that the minimal cost of reducing carbon emission is no more than the cost of doing nothing at all.

Uniqueness. Define
Comparing the definition with that of M 1 in subsection 3.2, we find that it is needed for solutions to decrease exponentially to zero as x → −∞ such that the following uniqueness theorem holds: The risk-free interest rate is the current annual RMB deposit rate in China: r = 0.0325 3 .

5.2.
Marginal carbon abatement cost curve. Marginal abatement cost curves (MACCs) are used to obtain the carbon abatement cost function, g 1 (C) and g 2 (A), in this model. A marginal abatement cost (MAC) curve usually indicates the relationship between the marginal abatement cost (the cost to reduce the last unit of emission) and the total amount of emission reduction (Kesicki [29]; Morris et al [36]; Criqui et al. [13]). Since the marginal cost is the derivative of total reduction cost with respect to the reduction amount (Jakob [27]), the total cost of reducing a given amount of carbon emission can be obtained by integrating the MAC with regard to the amount (Kesicki [29]). There are MACCs whose input is the percentages relative to the emission baseline as well, see Ellerman and Decaux [16]; Klepper and Peterson [32].
There are two approaches to obtain an MACC: expert-based approach and modelderived approach (Kesicki [29]).
In this paper, we mainly focus on the economical impacts of carbon emission abatement and the macroscopic decisions. Thus we choose top-down MACs to calculate the carbon abatement costs. Although the estimation of MACs differ from model to model, the magnitudes of different top-down MACs are similar (Ellerman and Decaux [16]; Klepper and Peterson [32]; Criqui et al. [13]). We choose the data obtained by the EPPA model from [13] to fit the MAC function 4 , and adjust the data measured in 1990 US dollars in that article to 2014 US dollars by multiplying the inflation factor 5 . See Appendix A for the adjusted marginal abatement costs associated with the amount of reduction of China. The EPPA MAC curve for China in [13] is shown in Figure 2 (a).  2), f (x) is an increasing function of x. We try to use the following two models to fit the MAC curve with the above properties: The data is approximated from Fig. 7, where the curve named China-MIT represents the EPPA curve for China. Although [32] Table 2 offers some accurate data from simulation, [13] Fig.  7 provides more points, making it more suitable for fitting. The measure of abatement amount in Fig. 7 (million tons of C, or MtC for short) is converted to billion tons of CO 2 , or GtCO 2 for short. 5 Source of US annual inflation rate from 1990 to 2013: http://www.statista.com/statistics/ 191077/inflation-rate-in-the-usa-since-1990/. The first model is parabolic model, and the second one is similar to the CRED model in [1,2], where m 2 is the upper bound for the potential amount of abatement, and the marginal cost approaches infinity as the reduced amount approaches m 2 .
The results of the fittings are shown in Table 1. The R 2 (coefficient of determination) values in Table 1 and Figure 2 (b) imply that the parabolic model outperforms the CRED model. Then the total abatement cost function g 2 (A) can be obtained by integration.

Commitment and penalty.
China doesn't belong to either the Kyoto Protocol or the EU ETS coverage where the commitment targets are imposed. Meanwhile, the voluntary carbon reduction target announced by China concerns about CO 2 emission intensity per unit of GDP instead of total emission amount (Teng and Xu [42]). As a result, we have to choose the reduction target and non-compliance penalty of EU ETS to approximate the case with China in this article.
Until now, three operational periods of EU ETS have been implemented or agreed on: Phase I from 2005 to 2007, Phase II from 2008 to 2012, and Phase III from 2013 to 2020 (Carmona et al [6]; International Emissions Trading Association (IETA) [26]; European Commission [12]). The countries belonging to EU ETS have to reduce their carbon emissions by 20% compared to their 1990 levels by 2020(IETA [26]; European Commission [11]), which is 6 years from now. Then the emission threshold isĪ = 0.8 × 2.4894 = 1.9915(GtCO 2 ). The non-compliance penalties for Phases I and II are e40/ton and e100/ton respectively(Carmona et al [6]; Cetin and Verschuere [8]), and from 2013 the penalty will increase annually with the annual rate of inflation in the Eurozone(IETA [26]; Fallman [18]; European Commission [11,12]). For simplicity, we choose the penalty in Phase II, i. e. e100/ton (about $136.5/ton). Then the parameters in our model are determined as follows: a 0 = 0.1365(thousand dollar/ton).
6. Numerical simulation. The parameters used in this section are the same as is defined in (5.1,5.3).
6.1. Semi-implicit method. In this section, we use finite difference method to solve the HJB equations. There are two frameworks of finite differencing: the explicit scheme (Song [41]; Fleming and Soner [21]) and the implicit scheme (Wang and Forsyth [44]; Forsyth and Labahn [22]). The time step in an explicit scheme is limited to guarantee the stability of the scheme, while in an implicit scheme there's less limitation on the time step. Thus, the implicit scheme is usually faster. Forsyth ( [44,22]) proposed an implicit iteration method to solve HJB equations. However, the method in [22] didn't work well on the equations in this paper. Therefore, a semi-implicit method is applied.

Recall equation (3.3) in section 3. Let y = x − lnĪ, and equation (3.3) can be turned into
Consider the following linear problem on Q mT Ω m ×[0, T ], Ω m [−m lnĪ, m lnĪ]: At time t j , we use the optimal policy at time t j−1 and the implicit central/forward /backward differential scheme in [44,22] to calculate Φ j . Then C j can be obtained from the following definition: Then C j can be used to calculate Φ j+1 . Since the solution Φ to equation (3.3) belongs to C 2+α,1+α/2 (Q ∞,T ), C j−1 is a good approximation of C j with ∆t small enough. The process of numerically solving equation (4.1) is similar. We call this scheme semi-implicit because Φ j is calculated implicitly, and C j−1 , instead of C j , is used to calculate Φ j . From Figure 3, the scheme has first order convergence rate.  Figure 4. Relationship between carbon emission, time and optimal policy 6.2. Optimal policy analysis. In this section, we will present figures of the optimal policy A * in model 2, and discuss the properties of A * shown in these figures. Figure 4 implies that the optimal policy increases with initial emission amount y = ln(I/I 0 ), and reaches the upper bound when y is large enough.
The relationships between the optimal policy A * and parameters are shown in Figure 5. These figures imply that the A * increases with µ, σ and ρ, and decreases with r.
From the parameters obtained by linear regression, we can infer that the growth rate of GDP (µ) and population (ρ) is positively correlated to the growth rate of carbon emission I (a 1 > 0, a 2 > 0). As a result, higher values of µ or ρ leads to higher CO 2 emission growth rate. This increases the probability of CO 2 emission exceeding the threshold at time T . In order to avoid the high penalty, policy makers decide to take a stricter reduction policy, as is implied in Figure 5 (a)(d). As σ increases, the volatility in CO 2 emission growth rate increases, and the uncertainty in carbon emission is higher. In order to make sure that the emission at time T won't be too high, the decision makers have to choose a more aggressive policy. Meanwhile, the relationship between C * and interest rate r captures the timing strategy in decision making. Higher interest rate r means a higher future value of the present cost of carbon abatement, or equivalently, a lower present value of future carbon abatement costs. From the point of reducing total costs in the whole process, the policy makers tend to lower the current policy, or delay a strict policy. This explains Figure 5 (c).
From Figure 4 it can be inferred that the space R × [0, T ] can be divided into two regions, depending on whether the upper bound of control policy is reached: and Σ 2 lies on the right side of Σ 1 . We use the policy boundary, i.e. the boundary between Σ 1 and Σ 2 , to capture the impact of the upper bound of the control policy,   It can be inferred from Figure 6 (a) that, although the policies near B(t) differ with different upper bounds, they are the same when the values of optimal policies are relatively small, and the policies will converge to the one withoutĀ asĀ grows larger. Figure 6 (b) implies that the boundary B(t) decreases with t, and increases withĀ. AsĀ grows, B(t) moves to the right side (resulting in larger boundary values) gradually, and approaches infinity asĀ goes to infinity.
The properties of the optimal policy in Model 1 are similar.
6.3. Minimal cost analysis. In this section, figures of the minimal cost changing with different parameters will be shown, and properties of the minimal cost will be analyzed. The effect of the upper bound on the minimal cost is also discussed. Figure 7 shows the minimal cost Ψ changing with emission y = ln(I/Ī) and time t. It can be inferred that Ψ is an increasing and concave function of y. Figure 8 presents the relationships between the minimal cost Ψ and parameters. It can be inferred that Ψ increases with µ, σ and ρ, and decreases with r.  The reason for Ψ increasing with µ and ρ is that they both contribute to the carbon emission growth rate. At the same time, the country chooses a more aggressive policy to reduce carbon emission when µ and ρ is higher, resulting in a higher abatement cost, see Figure 5 (a)(d). The case with r is on the contrary. Higher values of r leads to a smaller discounted value of future costs as well as a loose control policy, which make the total cost lower. Since the optimal policy A * increases with the volatility σ, the total cost increases as well.
Since both GDP and population growth rates are positively correlated with carbon emission growth rate, and high volatilities in GDP growth lead to higher emission reduction total costs, a country should prevent its population and economy from developing too fast, and keep a relatively steady economical growth, if it want to reduce the emission reduction cost. For given GDP and population growth rates,  higher a 1 or a 2 (representing the contributions of these growth rates in increasing the emission growth rate) leads to a higher total cost. If decision makers can take new techniques to reduce these values, they can reduce the total carbon abatement cost without slowing down economical growth. Thus technological improvement is crucial in this process. Figure 9 shows different minimal costs with different upper bounds for control policy. It is shown that the minimal cost Ψ decreases as the upper boundĀ increases, and Ψ converges to the minimal cost function without the upper bound as A goes to infinity.
The properties of the minimal cost in model 1 are similar.
6.4. Technical influence. Although the influence of technology isn't included in our model directly, it is reflected by some parameters. The parameters a 1 and a 2 in equation (2.3) represent the contribution of GDP and population growth rate in increasing CO 2 emission growth rate. m 1 and m 2 in the parabolic MAC fitting model (5.2) affect the costs to reduce carbon emissions directly. It is reasonable to believe that these values will be lowered if more improvements are made in clean development techniques, and the impacts of these parameters are a good representation of those of technology.  Figure 10. Relationship between optimal policy and technical parameters Figure 10 shows that the optimal policy increases as a 1 or a 2 increases, and decreases as the marginal abatement cost parameters, m 1 or m 2 increases. This implies that a country tends to take stricter control policies with high carbon emission growth rate, and it prefers more conservative abatement policies if the cost to reduce carbon emission is high. Meanwhile, Figure 11 shows that the total carbon reduction cost increases if a 1 , a 2 , m 1 or m 2 increases. This means that high emission growth rates lead to higher total carbon reduction costs. Meanwhile, when the marginal abatement costs are high, the total cost to reduce carbon emission will go up, disregard of the conservative control policy.  Figure 11. Relationship between minimal cost and technical parameters 7. Conclusion. In this paper, we analyzed the carbon emission reduction cost minimization problem of a country. We have established models of reducing the carbon emission growth rate as well as reducing the emission amount, which lead to the associated HJB equations. The existence and uniqueness of the solutions are established by PDE techniques. Then we used MAC curves and data to conduct numerical simulations, and discussed the properties of the total minimal costs and the optimal control policies.
This paper established a model with a guideline for a country to reduce its carbon emission with a relatively low cost. It also helps policy makers to make decisions of carbon abatement. The results in this paper imply that high GDP or population growth rate or high volatility in GDP growth leads to more aggressive carbon abatement policy and higher total cost. Technical parameters play an important role in affecting the total cost as well. For a country aiming at reducing carbon abatement costs in the long run, it is better to keep a relatively steady economical growth, control the population growth, and seek for more environmentally-friendly technology.