Evaluating the Causal Effects of Emissions Trading Policy on Emission Reductions Based on Nonlinear Difference-In-Difference Model

Abstract

Based on panel data from 30 provinces, cities, and autonomous regions from 2001 to 2019, this paper uses the nonlinear difference-in-difference ($DID$) method to estimate the distribution of causal effects of emissions trading policy on emission reduction in Chinese industrial enterprises, and examines the heterogeneity of the effects. The empirical results show that (1) the emissions trading policy has a significant effect on industrial $S{O}_2$ emissions reduction in China, where the reduction effect is larger in non-pilot areas than in pilot areas; (2) the policy effects are not proportional to the regional $S{O}_2$ emissions intensity, and the emissions trading policy is not more effective in regions with higher industrial $S{O}_2$ emissions intensities. One advantage of this paper is the use of nonlinear $DID$ to estimate the emissions reduction effect, which eliminates the bias problem caused by the strict linearity assumption of the classical $DID$ method. Another advantage is that the combination of the random forest method avoids the subjectivity in the selection of control variables and uses distribution effects for multilevel comparisons. This method improves the validity of estimating the effect of emissions trading policy and provides targeted policy suggestions for the effective promotion of system implementation, all of which have academic and application value.

1. Introduction

The rapid growth of China’s GDP since 1978 has made it the second-fastest-growing economy in the world in recent times [1]. However, environmental problems have become increasingly serious, with 135 Chinese cities exceeding air-quality and pollution standards and 34.0% experiencing acid rain in 2020, an increase of 0.7 percentage points from 2019 [2]. Preventing and controlling pollution and protecting the economic environment have become urgent tasks. The emissions trading system is the first large-scale market-based environmental regulation policy in China. It plays a regulatory role in resource conservation, environmental protection, and green development for environmental resource allocation, mainly using market mechanisms. Emissions rights trading originated in the United States, and the American economist Coase formally proposed pollutant emissions rights in 1960. China’s emissions rights began to be explored early on, but progress has been slow [3]. Since 1987, emissions trading has been carried out in several pilot cities, one after another. However, problems with the secondary market still persist for local emissions trading and when setting benchmark prices for emissions trading in China [4]. The effect of emission reduction trading on industrial $S{O}_2$ reduction has long been debated, so this issue in China is a fascinating and crucial subject to study.

China launched the “4+3+1" program in 2002, but due to the early implementation and lack of practical experience at that time, many pilot regions had zero trading volume. The pilot program was further expanded in 2007, and the pilot system was improved. Each pilot region introduced its own policies to suit the local conditions, and the scope and scale of the transactions also expanded [5]. Therefore, this paper constructs panel data for 30 provinces from 2001 to 2019 based on the 2007 pilot policy as a quasi-natural experiment. We use 11 approved provinces, such as Jiangsu and Tianjin, as the treatment group and 19 other provinces as the control group and adopt nonlinear $DID$ to evaluate the following questions: Has China’s first market-based environmental regulation policy had a positive effect on its emissions reduction? What is the extent of its emissions reduction? Are there differences in the policy treatment effects for the treatment and control groups? Is the implementation of the policy more effective for enterprises or regions with higher $S{O}_2$ emissions intensities?

2. Literature Review

Much research has been conducted on the emissions trading policy. Montgomery (1971), Woerdman (2000) demonstrated the policy’s feasibility and studied the relevant policy implementation systems, indicating that the policies should be implemented based on environmental conditions. Pollution discharge fees should be increased for systems with higher levels of pollution [6,7]. Cason (2003) proposed that the strict management of emission enterprises is insufficient to promote emissions trading effectively. At the same time, it is necessary to strengthen government oversight and emissions management [8]. Rausch and Abrell (2017) investigated a hybrid emissions trading system based on the uncertain corporate emissions reduction costs and future emissions. The findings show that a hybrid policy that introduces price or reduction boundaries allows for hedging differences in marginal reduction costs across subregions [9].

With regard to the emissions reduction effect of the emissions trading policy, many studies have been conducted in the related literature. Most scholars believe that the emissions trading policy is effective in reducing emissions. Yan et al. (2012), Li et al. (2016), and Zhang et al. (2017) conducted quasi-natural experiments with pilot implementations of the 2002 emissions trading policy in six provinces and cities of China. Their analyses discovered that the emissions trading policy is able to improve pollution reductions [10,11,12]. Using the $DEA$ model, Tu et al. (2018) discovered that the emissions trading system not only reduces pollutant emissions but also aids in economic development [13]. Tang et al. (2017) discovered that the emissions trading policy helps to reduce $S{O}_2$ emissions while reducing corporate profits [14]. Wu et al. (2018) discovered a weak “Porter effect” from emissions trading in China [15].

The DID method is the most commonly used method in this field. For example, Fu et al. (2018) and Wu et al. (2021) used the $DID$ method and the $DID-PSM$ method to find that the efficiency of green development can be improved by the $S{O}_2$ emissions trading policy, which has little effect on promoting green development [16,17]. Applying the $DID$ method, Shi et al. (2020) discovered that the system’s implementation positively impacts energy consumption and improves the green total factor energy efficiency [18]. Ren et al. (2019) and Fu et al. (2019) used the $DID$ method to evaluate the effects of further increasing the policy pilots after 2007. They found that the implementation of the policy has shown positive effects for both $S{O}_2$ emission reduction and economic growth, which is a win–win in terms of environmental and economic goals [19,20]. Si et al. (2020) found significant emission reduction effects of emissions trading policy on pollutants $S{O}_2$, $N{H}_3$, and $COD$ [21]. Qi et al. (2020) evaluated the effects of emissions trading policies by using the $PSM-DID$ method. The results showed that there was a significant emission reduction effect of policy on industrial $S{O}_2$ and industrial wastewater emission results, but there was no significant contribution to green development in the short term ([22]. Yu et al. (2021), Cai et al. (2022), and Luo et al. (2018) conducted studies on carbon trading rights, demonstrating that carbon trading in China exhibits upward and then downward trends for the development of carbon technical efficiency [23,24,25]. With the use of Support Vector Regression (SVR), Linear Regression (LR), and Analytic Hierarchy Process (AHP), Hong et al. (2022) conducted data mining and the evaluation of pollutant emissions, market prices, marginal benefits of emissions rights, and initial emissions rights [26]. Qi et al. (2021) implemented an emissions trading policy using multi-period $DID$, which indicates the policy significantly increased Chinese enterprises’ OFDI [27]. Using the $DID$ model, Wu et al. (2022) discovered that actions towards the prevention and control of air pollution can significantly improve air quality in vital resource-based cities [28]. Zhang et al. (2020) used the $DID$ model to examine the impact of the emissions trading system on environmental efficiency. The results indicate that the policy has significantly improved the environmental efficiency in these pilot provinces [29]. Li et al. (2022) explored the impact of an emissions trading system (ETS) on the regional industrial structure. Their results show that the implementation of ETS can promote upgrades in the regional industrial structure but can hinder the rationalization of these upgrades [30].

Another body of literature contends that emissions did not decrease with the implementation of the emissionstrading system. Zhang et al. (2020) used the $DID$ method to study the environmental and employment impacts of emissions trading policy. They found that piloting the policy did not improve the environment, but significantly increased employment levels [31].

As mentioned by Xiang et al. (2016), there is a growing body of literature on natural experiments on policy implementation using the $DID$ method due to its simplicity and ease of implementation. However, this method also has many shortcomings. One of the biggest criticisms is that it only elaborates on linear relations well and no longer works for nonlinear relations. When using the $DID$ method, some studies have obtained biased results because of ignoring the strict assumptions required [32]. Another shortcoming is that the method can only give average treatment–effect estimates and not other quantile estimates. Athey and Imbens (2006) developed a nonlinear $DID$ method, referred to as “changes-in-changes” ($CIC$), to assess the impact of policy changes that replaces the widely used $DID$ model. Unlike traditional methods, the authors proposed a nonlinear model that allows unobservable effects to change over time. This nonlinear $DID$ method estimates the entire counterfactual distribution of outcomes between the treated and untreated groups [33]. Lucas and Mbiti (2008) investigated the relationship between Kenya’s compulsory primary education policy and entrance examination grades using the nonlinear $DID$ model. The inhibitory effect of student attendance was significant at the 1% level, with a positive but small effect at the middle quantile level and a positive effect at the higher quantile [34].

The literature on the effect of emissions trading policies on emission reductions has three main characteristics. First, the current research primarily employs the $DID$ method, but in actual natural experiments, meeting the strict requirements of the $DID$ method is challenging. However, the results are assumed to be accurate. Second, current literature assumes that the application of the policy will have the same treatment effect on various groups, even though various experiments should be conducted to account for individual differences. Third, while all previous studies have been able to determine the average treatment effect of the policy, they have not been able to determine the distribution treatment effect. However, the distribution treatment effect explains which geographic areas the policy is effective for. In this study, we construct panel data of 30 provinces during the period from 2001 to 2019. The 11 approved provinces of emissions trading policy are evaluated as the treatment group, and the other 19 as the control group using the nonlinear $DID$ model. The significance of this work is as follows. (1) We use the nonlinear $DID$ method to overcome the problem of errors caused by the fact that natural experiments do not satisfy the strict assumption of $DID$. (2) Due to differences in the characteristics of different groups, we expect a different effect. (3) The treatment effect is further obtained to analyze whether a more obvious emissions reduction effect occurs with the implementation of the policy in regions with higher $S{O}_2$ emissions intensities.

3. Materials and Methods

3.1. The Classical DID Method

For $i=1\dots N$, the data are presented as $(Y_i,G_i,T_i)$. We use the indicator $G_i\in \left\{0,1\right\}$ to denote the treatment group ($G_i=1$) and the control group ($G_i=0$), where $T_i\in \left\{0,1\right\}$ is a treatment variable with $T_i=1$ if individual i receives treatment, and $T_i=0$ if individual i does not receive treatment. We use $Y_i$ to denote the outcome and $Y_i^I$ and $Y_i^N$ to denote the outcome of individuals with and without treatment, respectively. Let $I_i=G_i\times T_i$; we have

$$\begin{array}{c}\hfill Y_i=Y_i^N\times (1-I_i)+I_i\times Y_i^I.\end{array}$$

(1)

In the classical $DID$ model, the outcome of individual i without the intervention satisfies the following:

$$\begin{array}{c}\hfill Y_i^N=\alpha +\beta \times T_i+\eta \times G_i+{\epsilon }_i,\end{array}$$

(2)

where $\beta$ and $\eta$ represent the time effect and group effect, respectively, and ${\epsilon }_i$ is a random error with ${\epsilon }_i\perp (T_i,G_i)$. Classical $DID$ estimates are obtained using the following:

$$\begin{array}{c}\hfill {\tau }^{DID}=E\left[Y_i\mid G_i=1,T_i=1\right]-E\left[Y_i\mid G_i=1,T_i=0\right]\\ \hfill -[E\left[Y_i|G_i=0,T_i=1\right]-E\left[Y_i|G_i=0,T_i=0\right]].\end{array}$$

(3)

3.2. The Nonlinear Difference-In-Difference Method

The $DID$ method usually focuses on the average causal effect of the treatment group. As the policy effect varied with unobservable features and the distribution of individuals in the group are different, the average treatment effect on the two groups is not the same. Although the classical $DID$ method requires relatively few assumptions to calculate the effect of policy interventions on the treatment group, the effect of policy interventions on the control group is rarely considered, which is a disadvantage of the classical $DID$ method. Therefore, there is controversy in the literature about the results of using classical $DID$ method. Athey and Imbens (2006) proposed the changes-in-changes ($CIC$) model—a nonlinear $DID$ method—for estimating the effects of policy interventions on control groups. In particular, the method assumed that a treatment effect depends on the unobservable characteristics of an individual rather than directly on the population. The distribution of treatment effects differed due to the different distribution characteristics of the groups [33].

3.2.1. The Changes-In-Changes Model

First, a random variable $U_i$ is introduced to represent the region’s unobservable characteristics in order to better identify the $CIC$ model. To simplify the analysis, we drop the the subscript i and use $\left(Y,G,T,U\right)$ to denote the vector of random variables. Then, we define the following:

$$Y_{gt}^N=Y^N\mid G=g,T=t,Y_{gt}^I=Y^I\mid G=g,T=t,$$

$$Y_{gt}=Y\mid G=g,T=t,U_g=U\mid G=g,$$

where $Y_i=Y_i^N\times (1-I_i)+I_i\times Y_i^I$, and $I_i=G_i\times T_i$ represents the treatment indicator. $F_{Y^N,tg},F_{Y^I,tg},F_{Y,tg}$ and $F_{U,g}$ denote the corresponding distribution functions. $Y_{gt}$ represents the untreated subgroups, which are subgroups other than $(g,t)=(1,1)$. At the same time, a set of hypotheses about the distribution of the counterfactual status of the treatment groups in terms of the second-period outcomes are analyzed. Therefore, we can use the simultaneous distribution of the observable values $(Y,G,T)$ to express the distribution of $F_{Y^N,11}$. Indeed, these findings suggest that $F_{Y^N,11}$ can be expressed by conditional outcome distributions about the other three subpopulations $F_{Y,00},F_{Y,01}$, and $F_{Y,10}$, which are observable.

Assumption 1 (Model): In the absence of an intervention, the outcome of an individual satisfies the relationship $Y^N=h(U,T)$, where h is an unidentified, nonlinear function. Since h does not change as the group changes, differences between groups result from the distribution of U. To identify $F_{Y^N}$ with this model, the following hypotheses are required.

Assumption 2 (Monotonicity): For $t=0,1$, the production function $h(U,T)$ strictly increases in U.

Assumption 3 (Time Invariance): The structure of a group remains constant over time. Given a group G, U as an unobservable variable is time-stationary; that is, $U\perp T\mid G$.

Assumption 4 (Support): $4supp[U\mid G=1]\subseteq supp[U\mid G=0]$.

Different from the $CIC$ model, the classical $DID$ model has two additional assumptions:

$$U=\eta \times G+\epsilon ,h\left(U,T\right)=\alpha +\beta \times T+U,$$

that is, $Y_i^N=\alpha +\beta \times T_i+\gamma \times G_i+{\epsilon }_i$. As a result, the classical $DID$ model is a special case of the $CIC$ model, which relaxes and generalizes the classical $DID$ model’s assumptions.

Assumptions 1–3 are collectively referred to as the $CIC$ model, and the authors propose that Assumption 4 be invoked selectively as needed. According to Assumption 1, the results are not obtained directly depending on the group, and all associated unobservable variables must be encapsulated in a simple indicator U. According to Assumption 2, higher unobservable values must match higher outcomes. Weak monotonicity is only a standardization in a given subpopulation. Since we assume that unobservable values cause higher values in both periods, it is only limitative. When the unobservable variable is a personal characteristics such as health or ability, this structure naturally arises. Additive models allow for the existence of rich non-additive structures and can satisfy strict monotonicity automatically. As the outcome $Y_{gt}$ is continuous in models, the difference between strict and weak monotonicity is minimal. Furthermore, this assumption is overly restrictive if particles exist in the $Y_{gt}^N$ distribution. Assumption 3 requires that the population of a given group remains constant over time. The $CIC$ and $DID$ methods are based on the assumption that the differences between groups must be stable to achieve the state of affairs in which one group can be used to eliminate trends in the other. Assumption 4 implies that $supp\left[Y_{10}\right]\subseteq supp\left[Y_{00}\right]$ and $supp\left[Y_{11}^N\right]\subseteq supp\left[Y_{01}\right]$, and this assumption is relaxed in the authors’ subsequent inferences.

3.2.2. Treatment Effect on Treatment Group

Based on the above assumptions, the counterfactual distribution can be obtained after the following derivation:

$$\begin{array}{c}\hfill F_{Y^N,11}=F_{Y,10}\left(F_{Y,00}^{-1}\left(F_{Y,01}\left(y\right)\right)\right),\end{array}$$

(4)

where the distribution function of $F_{Y^N,11}$ is unobservable but we can derive the other three distribution functions from the data. Using these three observable distribution functions, we can obtain the unobservable distribution function. The average treatment effect $\left(ATT\right)$ can be written as follows after a series of transformations:

$$\begin{array}{c}\hfill {\tau }^{CIC}=E\left[Y_{11}^I-Y_{11}^N\right]=E\left[Y_{11}^I\right]-E\left[K^{CIC}\left(Y_{10}\right)\right]=E\left[Y_{11}^I\right]-E[F_{Y,01}^{-1}\left(F_{Y,00}\left(Y_{10}\right)\right].\end{array}$$

(5)

Clearly, no assumptions are made about the specific form of the h function in the preceding process, so the entire identification process is non-parametric. Moreover, the effect can be estimated using the sample mean and actual distribution. In other words, the $CIC$ model does not require assumptions as strict as those of the $DID$ model and generalizes the model.

3.2.3. Treatment Effects on the Control Group

To date, we have constructed a model with no intervention. In the presence of an intervention, no outcome model is required to infer the effect of the policy change on the treatment group, that is, the effect of the "treatment" on the treated. We simply need to compare the treatment group’s actual outcome to the counterfactual. However, more structure is needed to analyze the effect of policy intervention on the control group.

First, note that the difference between the counterfactual distribution of $Y_{01}^I$ and $Y_{11}^N$ is qualitative. After all, three subgroups do not receive treatment, all of which could be used in determining the outcome distribution without treatment in a fourth subgroup. Though only one subgroup received the treatment, we still want to know the distribution of $Y_{01}^I$. Therefore, we can build a transformation based on group 1 and apply it to $Y_{00}$, that is, assuming within a group that the distribution of U does not change over time. More specifically, $Y_{01}^I=h^I\left(U_0,1\right)$, $Y_{00}=h\left(U_0,0\right)$, and

$$\begin{array}{c}\hfill Y_{01}^I\sim h^I\left(h^{-1}\left(Y_{00};0\right),1\right).\end{array}$$

(6)

Because the distribution of $U_1$ remains constant over time, for $y\in supp\left[Y_{10}\right]$,

$$\begin{array}{c}\hfill F_{Y^I,11}^{-1}\left(F_{Y,10}\left(y\right)\right)=h^I\left(h^{-1}\left(y;0\right),1\right).\end{array}$$

(7)

The roles of group 0 and group 1 in $K^{CIC}\left(y\right)$ are simply switched. Following this logic, we apply the method described previously to calculate the counterfactual distribution of $Y_{01}^I$. In other words, replacing G with $1-G$ and taking Assumptions 1–3, we assume that $Y^I=h^I(U,T)$, where $h^I\left(u,t\right)$ increases strictly in U. Then, we can determine the distribution of $Y_{01}^I$ from the distributions of $Y_{00}$, $Y_{10}$, and $Y_{11}^I$ on the restricted branch set $\left[Y_{11}^I\right]$ as follows:

$$\begin{array}{c}\hfill F_{Y^I,01}\left(y\right)=F_{Y,00}\left(F_{Y,10}^{-1}\left(F_{Y^I,11}\left(y\right)\right)\right).\end{array}$$

(8)

3.3. Model Selection

To study the effect of emissions trading on industrial $S{O}_2$ emissions reduction, in this paper, we use the $CIC$ model. To be more specific, we begin with $h(U,0)=U$. In this case, U represents the industrial $S{O}_2$ emissions intensity, which is demonstrated by the regional emissions trading policy at period 0 after accounting for other control variables. Because of the impact of heterogeneity between groups and the effects caused by the emissions trading policy at period 0, the distribution of $U\mid G=g$ should have a discrepancy between different groups. In addition to the assumption that U is normalized in this case, the application of the $CIC$ model must satisfy two assumptions. First, the distribution of U among the group does not change with time. Because U represents the region’s characteristics, changes in the emissions trading policy cannot cause regions to adjust quickly and should not influence whether to implement the emissions trading policy. This assumption is reasonable. Second, in the absence of policy interventions, the resultant functions $h(U,1)$ of the two groups should be the same. This assumption excludes the relationship between the industrial $S{O}_2$ emissions intensity and the pilot area’s emissions trading policy over time. On this basis, two additional assumptions are required by the more stringent $DID$ model. (1) The main difference between the pilot and non-pilot regions of the emissions trading system is the industrial $S{O}_2$ emissions intensity. (2) The additive effect is the same for all individuals over the process of change over time. Based on natural quasi-experiments, these assumptions are difficult to prove. The distribution of industrial $S{O}_2$ exhibits different shapes with different groups and over different time periods, so we use the $CIC$ model here to investigate the effect of the emissions trading policy.

When applied to the problem studied in this paper, the treatment is the emissions trading policy. The outcome variable Y represents the $S{O}_2$ emissions intensity; G indicates whether the provinces, cities, or regions are pilot areas for emissions trading, presented as either 1 or 0; and T represents the time indicator variable of whether the pilot policy of emissions trading begins to be implemented, and is 1 if the start date is 2007 or later and 0 otherwise.

3.4. Variable Selection and Data Description

3.4.1. Explained Variable

The industrial $S{O}_2$ emissions intensity (lnpso2) is expressed as industrial $S{O}_2$/gross industrial output (tonnes/billion CNY), providing a more accurate description of industrial $S{O}_2$ emissions.

3.4.2. Core Explanatory Variables

Time, the dummy variables for each group, and the interaction terms of dummy variables for the policy implementation provinces and the policy implementation time are the main explanatory variables. Treat denotes the dummy variable for policy implementation provinces. If the province is one of the 11 provinces approved by the Ministry of Environmental Protection in 2007, the dummy variable Treat takes a value of 1; otherwise, it takes a value of 0. For the time dummy variable (Time), a start date of 2007 or later takes a value of 1, and that before 2007 takes a value of 0. The primary variable of interest is the interaction terms between the two dummy variables: Treat*Time.

3.4.3. Control Variables

The factors influencing the industrial $S{O}_2$ emissions intensity are controlled as thoroughly as possible based on data availability. The following 18 control variables are chosen: population size, industrial pollution control, investment in industrial waste gas control, sewage charges, economic development level, wage level, regional technological innovation level, technology introduction, technology level, fixed assets, industrialization level, industrial structure, scale of the service industry, degree of marketization, number of employees, scale of energy consumption, education development, and fiscal decentralization. At the same time, to carry out the relevant research in this paper, the industrial $S{O}_2$ emissions intensity, economic development level, wage level, regional technological innovation level, technology introduction, technology level, fixed assets, population size, investment in industrial pollution control, investment in industrial waste gas control, sewage charges, scale of energy consumption, and number of employees are all logarithmically treated based on the actual data.

By summarizing all the variables used in this paper,

Table 1. Definition and description of each variable and data source.

Variable	Symbol	Variable Description	Data Source	Units
Area Number	id	1–30	-	-
Time variables	year	2001–2019	-	-
Policy dummy variables	treat	0 or 1	Provincial and Municipal Environmental Protection Departments	Dummy variables
Time dummy variable	time	0 or 1	Environmental Protection Department of each province and city	Dummy variables
Industrial $S{O}_2$ emission intensity	lnpso 2	Industrial $S{O}_2$ emissions/Total industrial output	China Environmental Statistics Yearbook	Tonnes/Billion CNY
Regional technological innovation level	patent	Number of patent applications granted	China Statistical Yearbook	Individual
Economic development level	gdp	GDP	China Statistical Yearbook	Million/Person
Wage level	wage	Average wages of employees	China Statistical Yearbook	CNY
Scale of the service industry	service	Value added tertiary sector/Gross Domestic Product	China Statistical Yearbook	Score
Technology introduction	fdi	Amount of foreign direct investment	China Statistical Yearbook	Billion CNY
Fixed assets	fassets	Investment in fixed assets	China Statistical Yearbook	Billion CNY
Educational development	edu	Average years of schooling in the workforce	China Statistical Yearbook	Year
Fiscal decentralization	fiscal	Fiscal expenditure per capita by province/ (Fiscal expenditure per capita in each province + Fiscal expenditure per capita in the central government)	China Statistical Yearbook	Score
Industrialization level	indus	Secondary sector value added/Gross Domestic Product	China Statistical Yearbook	Score
Scale of energy consumption	energy 1	Total energy consumption	China Energy Statistical Yearbook	Million tonnes
Energy consumption structure	energy 2	Coal energy consumption/Total energy consumption	China Energy Statistical Yearbook	Score
Population size	pop	Number of resident population at the end of the year	China Statistical Yearbook	Million people
Number of employees	labour	Number of labour force at the end of the year	China Statistical Yearbook	Million people
Investment in industrial pollution comtrol	invest 1	Investment completed in industrial pollution control	China Environmental Statistics Yearbook	Million CNY
Investment in industrial waste gas comtrol	invest 3	Investment in industrial waste gas treatment	China Environmental Statistics Yearbook	Million CNY
Sewage charges	charge	Amount of sewage charge levied	China Environmental Statistics Yearbook	Million CNY
Degree of marketization	market	Local government expenditure/Gross Domestic Product	China Statistical Yearbook	Score
Technology level	rd	Industrial Enterprises Above Scale Internal expenditure on science and technology	China Science and Technology Statistical Yearbook	Million CNY

illustrates each variable’s definition and the data source.

Descriptive analysis is performed on all variables, and

Table 2. Descriptive Statistical Analysis of Variables.

Variable	Obs	Mean	Std. Dev.	Min	Max
id	570	15.5	8.663044	1	30
year	570	2010	5.482036	2001	2019
treat	570	0.3666667	0.4823177	0	1
time	570	0.6315789	0.4828001	0	1
lnpso2	570	4.803911	1.328462	−1.4503	7.6214
edu	570	8.553047	1.187864	3.960645	12.68113
energy1	570	8.854377	0.7369401	5.69454	10.63079
energy2	570	0.4495641	0.1595189	0.0163233	0.883037
patent	570	9.076958	1.677473	4.248495	13.1757
pop	570	8.166417	0.7529524	6.259582	9.351927
indus	570	0.4719952	0.1726606	0.0835022	1.672889
service	570	0.5794321	0.4741081	0.0648192	2.595965
fassets	570	8.60031	1.256432	5.252709	10.98615
invest1	570	11.57812	1.128235	6.914135	14.16335
invest3	570	10.87706	1.325285	4.941642	14.06343
charge	570	10.41566	1.037038	6.763654	12.56843
market	570	0.2099161	0.0949413	0.0767083	0.6283544
fiscal	570	0.8077755	0.0893072	0.3417719	0.9949374
wage	570	10.38492	0.6950548	8.97563	12.307
labour	570	7.572334	0.8125459	5.631212	8.874903
rd	570	13.61397	1.477453	8.747828	16.95744
gdp	570	9.138368	1.125942	5.6984	11.5868
fdi	570	4.833506	1.733361	−1.1713	7.7457
lso2	570	9.300873	10.62118	−0.0037	80.9593

describes the statistical analysis of the data.

3.5. Random Forest to Select Control Variables

There are two methods for ranking the importance of control variables in random forests, mainly based on mean decrease impurity and mean decrease accuracy. Mean decrease impurity means that for each tree, the control variables are ranked according to their impurity, and finally, the whole forest is averaged. Accuracy reduction is achieved by reordering the values of a column of features and observing how much it reduces the accuracy of the model. It directly measures the effect of each feature on the prediction accuracy of the model. For significant control variables, this method reduces the accuracy of the model significantly, while for insignificant control variables, this method has little effect on the accuracy of the model. Correspondingly, the two representations of feature importance ranking in the regression tree are %IncMSE and IncNodePurity. %IncMSE is the increase in MSE. Specifically, each control variable such as $X_1$ is assigned a random value. If $X_1$ is important, the prediction error will increase. So the increase in error is equivalent to the decrease in accuracy. The same is true for IncNodePurity. In the regression problem, node purity is actually the reduction of RSS (residual sum of squares). Node purity increase is equivalent to the reduction of the Gini coefficient; that is, the data in the node or classification are the same, known as Mean Decrease Gini.

Because the dependent variable in this paper is a continuous variable, a random forest regressor is constructed using the R package random forest. The number of trees used after parameter tuning is 500, and the leaf node defaults to one-third of the control variables. Based on the already constructed random forest regression model, the 18 control variables that influence the industrial $S{O}_2$ emissions intensity can be ranked in order of importance. Then, the cross-validation method is used to select the top-ranking control variables, which have a greater influence on the industrial $S{O}_2$ emissions intensities.

Figure 1 depicts the importance ranking of the 18 control variables in this paper by using random forests. The term “%IncMSE” represents the increase in mean squared error. By substituting random variables at random and observing the error in the model after the substitution, the variable becomes significant as the error increases. That is, the significance of the control variable is proportional to its magnitude; “IncNodePurity” denotes an increase in node purity. The sum of squares of residuals is used to calculate this indicator. Similarly, the magnitude of this variable is positively correlated with its importance. We discovered a discrepancy in the order of importance of control variables based on the different criteria of “%IncMSE” and “IncNodePurity.” This paper chooses the first few control variables with “%IncMSE” as the standard based on the specific situation that affects the industrial $S{O}_2$ emissions intensity.

To improve the model’s performance, determining how many control variables are best suited for model prediction is critical. In this paper, the number of control variables to be added to the model is determined by ten-fold cross-validation. The smaller the cross-validation error, the better the model fits. At the same time, the number of control variables is optimal. Figure 2 depicts a ten-fold cross-validation graph, with the horizontal axis representing the number of feature variables and the vertical axis representing the cross-validation error. We discover that adding more irrelevant control variables does not improve prediction accuracy, so having more control variables is not better. The cross-validation curve above clearly shows the relationship between model error and the number of control variables used for fitting. The error increases when the number of control variables is 0–2, but it gradually decreases when the number of control variables is 2–5. It is probably smaller when the number of control variables is 5. As the number of control variables increases, the error becomes larger. As a result, we select five control variables to minimize error and achieve the ideal state. Therefore, using “%IncMSE” as the control variables sorting criterion, the top five essential control variables are chosen to add to the model, which are wage level, energy consumption structure, regional technological innovation level, foreign direct investment level, and education level.

4. Results

This paper evaluates the emission reduction effect of the emissions trading policy on industrial $S{O}_2$. The difference in the policy treatment effects between the treatment and control groups is also confirmed. Moreover, the study illustrates the effect of policy implementation on emissions reduction for enterprises or regions with higher $S{O}_2$ emissions intensities.

4.1. Empirical Analysis of the Effect of Industrial $S{O}_2$ Emissions Reduction

We use a random forest approach to select the control variables needed in the model, which are wage level, energy consumption structure, regional technological innovation level, foreign direct investment, and education level. Previous studies assumed that policy implementation had the same treatment effect on the treatment and control groups. However, due to individual differences, the treatment and control groups had different treatment effects. This paper primarily employs the $CIC$ model based on relevant panel data from 30 provinces and cities from 2001 to 2019 with the $CIC$ command in STATA16.0 to make an empirical analysis.

Table 3. Treatment effects.

Variable	Treatment Group	Contral Group
Treat*time	−0.1702 **	−0.2968 **
	(0.0918)	(0.1301)
wage	−0.6109 ***	−0.6109 ***
	(0.1464)	(0.1267)
energy2	1.6652 ***	1.6652 ***
	(0.1874)	(0.1931)
patent	−0.2869 ***	−0.2869 ***
	(0.0411)	(0.0392)
fdi	−0.0512 *	−0.0512 *
	(0.0252)	(0.0267)
edu	−0.1910 ***	−0.1910 ***
	(0.0347)	(0.0318)

Note: () represents the standard error term, *, **, *** represent significant at the 10%, 5%, and 1% levels, respectively.

shows the empirical results for the estimation of the policy effects, and the values in parentheses are the estimated standard errors obtained after 50 iterations of the self-help algorithm. With the highest importance as screened by random forest, these results are obtained by sequentially including the top five control variables, which are wage level, energy consumption structure, regional technology innovation level, foreign direct investment level, and education level.

From the coefficients of the control variables, the magnitude of the coefficient of wage level is −0.6109, which shows that an increase in wages reduces the emissions of industrial $S{O}_2$. Wage growth indicates an improvement in the enterprises’ development level. To pursue higher-quality development, enterprises will pay more attention to the scale of pollutants generated and minimize the pollution of the environment while gaining benefits. The state can establish relevant policies that aid in the development of enterprises and raise the development level of enterprises to reduce emissions. The coefficient of the energy consumption structure is 1.6652, which means that more pollutants will be released as it rises. This will have a negative effect on industrial $S{O}_2$ emission reduction. In order to limit pollutant emissions, coal should be burned as infrequently as possible and fossil fuel energy should be substituted with clean new energy. The regional technological innovation level coefficient is −0.2869. In this study, we analyze the number of granted patent applications as a proxy for the degree of innovation in the region. We find that regional technical innovation is positively correlated with the ability of industrial $S{O}_2$ emissions reduction. As regional technological innovation levels increase, businesses use new technologies to carry out production tasks and cleaner technologies to address pollutant emissions with the achievement of emissions reduction. Foreign direct investment has a coefficient of −0.0511. The level of technology introduction in this paper is represented by foreign direct investment. It can be seen in Table 3 that foreign direct investment can promote the effect of industrial $S{O}_2$ emissions reduction, but only marginally. Foreign investment can improve production technology and increase technological investment, thus improving production technology and reducing industrial $S{O}_2$ emissions. The coefficient of education level is −0.1910, indicating that increasing the average education level has an inhibitory effect on pollutant emissions. Improving citizens’ education levels positively affects people’s quality of life, which promotes human capital accumulation and reduces industrial $S{O}_2$ emissions.

4.2. Distribution Effect

In terms of the average treatment effect, the emissions trading policy shows a significant suppression effect on the industrial $S{O}_2$ emissions in both the pilot and non-pilot areas. Furthermore, the average treatment effect in the non-pilot areas is higher than that in the pilot areas. However, there are some subjective factors while selecting pilot areas. The higher the industrial $S{O}_2$ emissions intensity, the more likely the pilot areas are selected. Therefore, the 10th and 90th percentiles are added to the quartiles to analyze at which stage the emissions trading policy is more effective for industrial $S{O}_2$ emissions intensity.

As with the treatment in the previous subsection, the average treatment effect indicates that the emissions trading policy has a significant suppressive effect on the industrial $S{O}_2$ emissions intensity in the pilot areas. The distribution effect allows us to observe more precisely at which stage the industrial $S{O}_2$ emission intensity is positively affected by the policy implementation. As can be seen from

Table 4. Distributed treatment effects.

Subsites	Treatment Group	Contral Group
Mean	−0.1702 **	−0.2968 **
	(0.0918)	(0.1301)
$q10$	−0.6638 **	−0.9136 ***
	(0.2619)	(0.2460)
$q25$	−0.2163	−1.1279 ***
	(0.1357)	(0.4205)
$q50$	−0.0714	−0.1756
	(0.1128)	(0.1552)
$q75$	0.1043	0.0859
	(0.0795)	(0.0916)
$q90$	0.0199	0.0199
	(0.1048)	(0.0863)

Note: () represents the standard error term, **, *** represent significant at the 5%, and 1% levels, respectively.

, the implementation of the emissions trading policy has a significant suppressive effect on the industrial $S{O}_2$ emissions intensity at 10 percentile. The absolute values of the contral group coefficients are larger than those of the treatment group, indicating that the implementation of the policy has shown a stronger effect on emissions reduction in non-pilot areas. Only for industrial $S{O}_2$ emission intensity at the 25th percentile, the treatment effect for the control group is negative at the 1% level, with a coefficient size of −1.1279. The coefficients for industrial $S{O}_2$ emissions intensities at the higher percentile do not pass the significance test.

When the control variables are added sequentially, the results show that the policy is more effective at the lower percentile and that none of the coefficients at the higher percentile pass the significance test. This indicates the implementation of the emissions trading policy is less effective for “dirtier” industrial enterprises or regions but more effective for the regions with lower $S{O}_2$ emissions. However, these industrial enterprises or regions with higher pollutant emissions intensities are the main source of environmental damage. Therefore, the government should improve the market-trading mechanism in regions with lower industrial $S{O}_2$ emissions to further implement and promote the policy. Supervision and management should be strengthened in regions with higher industrial $S{O}_2$ emissions. In this way, the policy implementation can have significant effects on these enterprises, thereby achieving the goal of environmental protection.

5. Robustness Check

5.1. Placebo Test

For the robustness check, we replace the explained variable industrial $S{O}_2$ emissions with domestic $S{O}_2$ emissions. The robustness test results are shown in

Table 5. Placebo test.

	Treatment Group	Contral Group
Mean	−0.2180	0.2356
	(1.0017)	(1.2336)
wage	−2.9476 ***	−2.9476 ***
	(0.9029)	(0.8052)
energy2	35.5129 ***	35.5129 ***
	(4.2869)	(4.4650)
patent	1.9950 ***	1.9950 ***
	(0.4304)	(0.4142)
fdi	(0.4304)	(0.4142)
	(0.4358)	(0.3899)
edu	0.0930	0.0930
	(0.2928)	(0.3342)

Note: () represents the standard error term, *** represent significant at the 1% levels, respectively.

.

Since the emissions trading policy mainly targets industrial $S{O}_2$ emissions, no effect should be seen on domestic $S{O}_2$ emissions. According to the results, the average treatment effects of treatment group and control group are not significant, indicating that the implementaiton of the emissions trading policy has no significant effect on domestic $S{O}_2$ emissions regardless of the inclusion of the control variables. The results confirm the implementation of the emissions-trading policy plays a catalytic role in industrial $S{O}_2$ emissions reduction. Moreover, the estimation results are robust.

5.2. Replacement Policy Intervention Time Point

This paper uses 2007 as the policy-impact point, but some areas in China began to implement the emissions trading policy as early as 2002. As a result, 2002 is chosen as the start date of policy intervention to generate a new start date dummy variable, Time1 (if the start date is before the year 2002, the value is 0; otherwise, it is 1).

Table 6. Replacement policy intervention time point.

	Treatment Group	Contral Group
Mean	−0.3329	0.5423
	(3.3377)	(4.2634)
$q10$	−1.5500	2.8954
	(4.7020)	(10.5034)
$q25$	−2.7905	0.8571
	(3.0349)	(2.9364)
$q50$	−1.6975	−3.1183
	(2.5386)	(3.6534)
$q75$	−0.5105	11.7543 *
	(3.5240)	(4.3629)
$q90$	2.9460	5.4399 ***
	(11.8022)	(3.6034)

Note: () represents the standard error term, *, *** represent significant at the 10% and 1% levels, respectively.

reports the estimates of the emissions trading policy effects after replacing the start date of the policy intervention with 2002. From test results, the average treatment effect after replacing the start date dummy variable is clearly not significant in either the treatment or control groups. Furthermore, we see that the treatment effect of the emissions trading policy after replacing the start date dummy variable is also insignificant at the lower percentile. This indicates that the emission reduction effect produced by the emissions trading policy is not caused by the start date.

5.3. Excluding the Policy Impact of the Same Period

In response to the increasingly changing environmental issues, China has undertaken several pollution control policies since 2000. The more representative means of environmental regulation is the regulation on “the Administration of Emissions Fee Collection” adopted in 2002. The simultaneous implementation of multiple policies may lead to the estimation of the nonlinear $DID$ model confounding the effects of other policies. In this paper, we further exclude the effect of sewage charge collection on the estimation results by adding the total amount of sewage charges collected by provinces as a control variable. If the interaction with the pilot emissions trading policy is still significant after adding the control variable, the decrease in the scale of pollutant emissions is not entirely due to the collection of sewage charges.

Table 7. Treatment effect after excluding the impact of the pollution discharge fee collection policy.

	Treatment Group	Contral Group
Mean	−0.1889 ***	−0.2908 ***
	(0.0723)	(0.1127)
$q10$	−0.5082 ***	−0.9242 **
	(0.1609)	(0.3298)
$q25$	−0.2028 *	−0.8615 **
	(0.1182)	(0.3969)
$q50$	−0.0718	−0.0726
	(0.1052)	(0.1186)
$q75$	0.0503	0.0632
	(0.0826)	(0.0856)
$q90$	−0.0879	−0.1160
	(0.1254)	(0.1107)

Note: () represents the standard error term, *, **, *** represent significant at the 10%, 5%, and 1% levels, respectively.

reports the estimated effects of the emissions trading policy after removing the emissions charge levy policy. It can be seen from the table that the average treatment effects after adding the control variable of sewage charges on industrial $S{O}_2$ emissions intensity are −0.1889 and −0.2908, indicating that the effect of the emissions trading policy on industrial $S{O}_2$ emissions reduction is still significant. Furthermore, we look at its emissions reduction effect from the quantile points: the treatment effects for industrial $S{O}_2$ emissions intensity in the treatment group are −0.5082 and −0.2028 in the 10th and 25th percentiles, respectively, which are significant at the 10% level. In the control group, the treatment effects for industrial $S{O}_2$ emissions intensity are −0.9242 and −0.8615 in the 10th and 25th percentiles, respectively, which are significant at the 5% level. According to the treatment effect coefficients, the reduction effect of the control group is still greater than that of the treatment group. For the high quantile, the treatment effect was not significant in either the treatment or the control group. This demonstrates that the findings of this paper are sound.

6. Discussion and Conclusions

Using the panel data from 30 Chinese provinces and cities from 2000 to 2019, this paper employs the nonlinear DID method to analyze the of the emissions trading system’s effect on industrial $S{O}_2$ emissions reduction. Compared with the classical DID methods, there are two major advantages of the nonlinear DID method: first, it relaxes the linear assumption, which does not always hold in practice. Second, the method can give estimates of the distribution effects in different quartiles, whereas the classical DID method can only give an average treatment effect. In particular, this paper evaluates the distribution effects of industrial $S{O}_2$ emissions intensity at the 10th, 25th, 50th, 75th, and 90th percentiles. In addition, we verified the robustness by replacing the dependent variable with domestic $S{O}_2$ emissions, advancing the policy implementation time and excluding any contemporaneous policy effects. The main conclusions of this study are summarized as follows:

• First, a certain emissions reduction effect on Chinese industrial $S{O}_2$ emissions occurs with the implementation of the emissions trading policy. Increases in wage level, regional technological innovation level, foreign direct investment, and education level have significant inhibitory effects on industrial $S{O}_2$ emissions. The effect of the policy is slightly greater in non-pilot areas than in pilot areas.
• Second, the effect of the policy is slightly greater in non-pilot areas than in pilot areas. That is, the treatment effect of the control group is slightly larger than that of the treatment group. In terms of distribution effects, the areas with lower pollutant emissions intensity take the better policy’s emissions reduction effect. However, some subjective factors influence the selection of the pilot areas. The areas with higher pollutant emissions are more likely to be selected as pilot areas. As a result, the policy intervention has been demonstrated to not have the same treatment effect on every individual.
• Third, for both the pilot areas and the non-pilot areas, the implementation of an emissions trading policy has a significant emissions reduction effect on areas (enterprises) with low industrial $S{O}_2$ emissions intensities. However, areas with larger pollutant emissions did not show a significant effect. In fact, the emission reduction effect of the policy diminishes as the $S{O}_2$ emissions intensity increases. This also explains why the policy’s implementation has a slightly higher treatment effect on the control group than on the treatment group.

By reviewing the above conclusions, this study gives the following policy recommendations.

• The government should improve the market trading system in the pilot areas, unify the market trading standards, and strengthen the supervision of policy implementation in the pilot areas. The implementation of any policy requires clear laws and regulations as well as implementation standards. Therefore, the government should introduce and clarify the relevant regulations of the emissions trading market and unify the quota-allocation methods in each province and city. Moreover, the pricing and evaluation standards related to trading should be enhanced to promote the fair promotion of the emissions trading policy in the pilot areas.
• To ensure the smooth operation of emissions trading, the government should enforce relevant laws and regulations, strengthen the supervision of the emissions rights trading market in pilot areas, and impose appropriate penalties on illegal emissions units. In addition to the relevant regulations at the national level, the relevant regulations of provinces and municipalities should also be very clear to give full play to the positive role of the market in environmental management.
• The current emissions-trading policy for enterprises or regions with higher $S{O}_2$ emissions does not play a positive role in reducing emissions. Therefore, the government should strengthen the total volume control and quota allocation for areas with higher industrial $S{O}_2$ emissions and raise the trading fees. Total volume control and quota allocations are critical to the foundation, which are prerequisites for the emissions trading market’s smooth operation. For regions or enterprises with different levels of industrial $S{O}_2$ emissions intensities, a more reasonable total amount of pollutant emissions should be determined in combination with the local environmental capacity. Furthermore, a reasonable emissions reduction task should be prescribed according to the pollutant emissions intensity. At the same time, refining the list of enterprises and industries that conduct emissions trading and making a reasonable plan for the total amount of emissions in a targeted manner are suggested.