Does the “Blue Sky Defense War Policy” Paint the Sky Blue?—A Case Study of Beijing–Tianjin–Hebei Region, China

Improving air quality is an urgent task for the Beijing–Tianjin–Hebei (BTH) region in China. In 2018, utilizing 365 days’ daily concentration data of six air pollutants (including PM2.5, PM10, SO2, NO2, CO and O3) at 947 air quality grid monitoring points of 13 cities in the BTH region and controlling the meteorological factors, this paper takes the implementation of the Blue Sky Defense War (BSDW) policy as a quasi-natural experiment to examine the emission reduction effect of the policy in the BTH region by applying the difference-in-difference method. Results show that the policy leads to the significant reduction of the daily average concentration of PM2.5, PM10, SO2, O3 by −1.951 μg/m3, −3.872 μg/m3, −1.902 μg/m3, −7.882 μg/m3 and CO by −0.014 mg/m3, respectively. The results of the robustness test support the aforementioned conclusions. However, this paper finds that the concentration of NO2 increases significantly (1.865 μg/m3). In winter heating seasons, the concentration of SO2, CO and O3 decrease but PM2.5, PM10 and NO2 increase significantly. Besides, resource intensive cities, non-key environmental protection cities and cities in the north of the region have great potential for air pollutant emission reduction. Finally, policy suggestions are recommended; these include setting specific goals at the city level, incorporating more cities into the list of key environmental protection cities, refining the concrete indicators of domestic solid fuel, and encouraging and enforcing clean heating diffusion.


Introduction
Since the reform was launched in 1978, China has experienced rapid economic and social development; on average, over 9% of GDP growth has helped China become one of the biggest economies [1,2]. However, extensive developments have led to environmental predicaments [3], especially increasingly severe air pollution, which has aroused significant concerns from the public and the government [4,5]. In 2017, 70.137% of all 365 cities' ambient air quality exceeded the World Health Organization (WHO) Interim Target-1 (IT-1) of 35 µg/m 3 [6].
Air pollution is a major contributor to the global burden of disease [7], including cerebrovascular diseases [8], respiratory diseases [9], lung cancer [10] and mortality problems [11]. In China, about 1.2 million people died prematurely due to complication associated with air pollution in 2010 [12]. In 2019, it is estimated that 12% of all deaths were attributable to outdoor and household air pollution [13]. In view of the scale of the impacts of air pollution on public health [14], urgent action is needed to improve air quality. Air pollution has also become a central domain and a momentous bottleneck in realizing the sustainable development in China [15]. has made the BTH region the most serious air pollution region in China [3]. The average concentration of PM 2.5 in the BTH region was approximately 5.5-7.3 times higher than the safety standards (15 µg/m 3 ) of the WHO [31]; more worryingly, the situation would possibly be worse in winter and spring [32].
Many studies have pointed out that the concern and participation of the government (government documents and policies) is the most direct factor affecting environmental quality [33]. Therefore, it is essential to evaluate the impact of these different policies on the environment [34]. Based on the fact that cities should be central subjects of air pollution control, this paper aims to use high-resolution air pollution data to explore whether the effect of the BSDW policy varies in heating seasons, multiple forms of cities and cities in BTH region, so as to provide clear and constructive suggestions for the formulation of air pollution prevention policies in China in the future.
The structure of the article is as follows. In Section 2, we execute a brief literature review; Section 3 is the data description and methodology; Section 4 shows the results and discussion; Section 5 is the conclusion and policy implications.

Literature Review
Air pollution in the BTH region has attracted the attention of many researchers. Some researchers have performed analyses on the scientific causes of air pollution [35,36] and a number of studies have been conducted on the impact of air pollution [37][38][39][40]. Scholars usually performed econometric models to directly compare the changes of air pollution and test the effects of air pollution control policies, including the difference-indifference (DID) method [41,42], regression discontinuity design (RDD) [43], propensity score matching (PSM) [44] and the synthetic control method (SCM) [45]. Some recent studies have focused on the impact of the "Action Plan for Air Pollution Prevention and Control" [46][47][48] and the "clean heating" policy [49][50][51]. Huang et al. [52] used hierarchical Bayesian models and daily data to estimate the health impacts of the "Action Plan for Air Pollution Prevention and Control" from 2013 to 2017 in 74 key cities in China, and found that annual average concentrations of PM 2.5 , SO 2 and CO decreased by 33.3%, 54.1%, 28.2%, while no significant changes were seen in NO x (9.7% reduction; 95% CI −23.0 to 42.4) or O 3 (20.4% increase; 95% CI −30.1 to 71.0). Guo et al. [53] used the varying coefficient model to analyze the impact of clean electricity utilization on yearly air pollution in the BTH region and identify the long-term relationship between heating and air pollution. The results showed that the coefficient between coal consumption and pollutant emissions was positive, while electricity consumption was negatively correlated with pollutant emissions, which indicated that "replacing coal with electricity" in BTH region had played an active role in air pollution control. There were also analyses showing solicitude for key policy directions of pollution mitigation in the BTH region [54]. Jiang et al. [43] utilized the RDD method and found that the BSDW policy reduced the monthly average concentration of PM 2.5 and PM 10 by 14.49 µg/m 3 and 23.43 µg/m 3 . Xu et al. [55] used the GAINS IV Asia (Greenhouse Gas and Air Pollution Interactions and Synergies) model to assess the potential for yearly air pollution abatement, air quality improvement and associated costs of the BSDW policy in the BTH region.
The above researches provide us with very valuable ideas, and there are still some areas that can be improved. Firstly, few studies have involved a systematic analysis of the impacts of BSDW policy on air quality based on daily scale and monitoring point data [43,55]. Due to the lack of high-resolution geographical data of prefecture-level cities, even those studies based on monthly or daily data failed to identify the policy effect accurately. Therefore, a comprehensive evaluation of the effectiveness of the policy is urgently needed. Secondly, few studies have explored the effect of the BSDW policy in the heating seasons [43], leading to the failure of identifying the policy effect on this special period that had been proven to contribute to serious air pollution. Thirdly, there has been little analysis focusing on different policy effects on different cities [55] and air pollutants. As one of the important topics in the research fields of environmental economics and environmental analyses, it is necessary to assess the pure effect of the implementation of the BSDW policy, which will help to realize the modernization of the government governance system and capabilities of developing countries.

Data
Cities in the treatment group and control group are shown in Figure A1, which are in red and pink, respectively. The treatment group contains Beijing, Tianjin, Shijiazhuang, Tangshan, Handan, Xingtai, Baoding, Cangzhou, Langfang and Hengshui, which have begun to take the BSDW policy since 5 July 2018, while the control group includes Zhangjiakou, Chengde and Qinhuangdao, which are not included in the BSDW policy to take pollution control actions. We use China's high-resolution air pollution reanalysis data [56] from 1 January 2018 to 31 December 2018 to evaluate the impact of the BSDW policy on urban air quality systematically and accurately. The concentrations of pollutants and meteorological data of monitoring points in each city are extracted by ArcMap 10.8 and Python. Figure A1 also depicts all the monitoring points selected in this paper, i.e., 947 monitoring points. The treatment group includes 570 monitoring points. In this paper, the daily average concentrations of various pollutants (PM 2.5 , PM 10 , SO 2 , NO 2 , CO, O 3 ) are taken as dependent variables. Moreover, we present the policy effect on the AQI. In order to eliminate the interference of meteorological factors on the results, we specifically control the daily average temperature (temp) [57,58], daily average wind speed (u: longitude wind speed; v: latitude wind speed) [59,60], daily average relative humidity (rh) [61,62] and daily average pressure (psfc) [63] of each monitoring point.
As shown in Tables 1 and A1, the overall pollution level in BTH region is high, for example, the daily average PM 2.5 concentration of treatment group is higher than the first level concentration limit standard of GB3095-2012 (35 µg/m 3 ) and the standard of WHO's Air Quality Guidelines (AQG) (15 µg/m 3 ). Similarly, the daily average concentration of PM 10 is also significantly higher than 50 µg/m 3 (the first level concentration limit standard of GB3095-2012) and 45 µg/m 3 (WHO (AQG 2021)). The daily average concentrations of SO 2 , CO and O 3 are lower than WHO (AQG 2021) standard (40 µg/m 3 ; 4 mg/m 3 ; 100 µg/m 3 ) and first level concentration limit standard of GB3095-2012 (50 µg/m 3 ; 4 mg/m 3 ; 100 µg/m 3 ). For NO 2 , the daily average concentration of the treatment group is higher than WHO (AQG 2021) standard (25 µg/m 3 ) but lower than first level concentration limit standard of GB3095-2012 (80 µg/m 3 ). The maximum and minimum values of all air pollutants concentrations imply that the BTH region's air pollution has obvious regional differences. Figure 1 (derived from Table 1) shows the concentration distribution of air pollutants in the BTH region before (from 1 January 2018 to 4 July 2018; the air pollution concentration histogram is located above the horizontal axis) and after (from 5 July 2018 to 31 December 2018; the air pollution concentration histogram is located below the horizontal axis) the implementation of the BSDW policy. From Table 1 and Figure 1, we can draw a preliminary conclusion that the pollutant concentration in the treatment group decreased more than that in the control group after the implementation of the policy.

Empirical Model
To assess the impact of China's BSDW policy on the local air quality, this paper utilizes the DID method, which is widely used in the estimation of the causal effect of a treatment by using longitudinal data from the treatment and control group [64][65][66]. The DID model can control the systematic differences between the treatment and control groups and isolate the changes in the outcomes over time between the samples that were and were not affected by the policy [67]. It can also remove the biases that could be the result of trends caused by other factors (i.e., missing variables). The DID model can be defined as follows: time = 0, before policy implementation 1, after policy implementation (3) where the dependent variable Air quality i,t refers to the monitoring station's daily concentration of PM 2.5 , PM 10 , SO 2 , NO 2 , CO, O 3 and other air pollutant used in the robustness analysis (AQI). treat i,t represents the monitoring station i at day t. It is equal to 1 if the monitoring station is located in the city where the policy is implemented, otherwise it is equal to 0. time i,t is a dummy variable that represents the BSDW policy and is equal to 1 since 5 July 2018, otherwise it is equal to 0. time i,t × treat i,t notes the interaction item of these two dummy variables. X i,t represents control variables related to the environmental pollution, namely the meteorological factors (including daily temperature, daily longitude wind speed, daily latitude wind speed, relative humidity and daily surface pressure). The individual-fixed effects α i and the day-fixed effects δ t are included to exclude the interference of other possible confounders. β i is the vector of the coefficients of each independent variable. β 0 represents the constant term. β 1 and β 2 represent the coefficient of treat i,t and time i,t , respectively. The coefficient of time i,t × treat i,t (β 3 ) captures the reduction effect of the policy on the treatment group [68], namely the concentration reduction amount of air pollutants, which can effectively avoid the endogenous problems caused by omitted variable bias and accurately reveal the net policy impact [69]. The basic principle of DID method is to build the difference-in-difference statistics reflecting the policy effect, that is, the coefficient of the interaction term β 3 , by comparing the differences between the control group and the treatment group before and after the implementation of the policy, so as to obtain the desired net effect of the policy [70]. A negative estimated coefficient (statistically significant) denotes a significant reduction in the pollutant concentration in the region because of the BSDW policy, indicating that the emission-control policy has played a positive and significant role. ε i,t denotes random disturbance term.

Unit Root Tests and Hausman Tests
We proceeded with the unit root diagnostic test; namely, the Levin, Lin and Chu test was used to examine the stationary processes of our variables [71]. The results suggested that all variables were stationary at level. Then, we applied Hausman test to test whether individual effects were independent from other explanatory variables, so as to determine the applicability of the fixed effect model and random effect model [72]. The original hypothesis of Hausman test is that individual effects are independent from other explanatory variables, namely the random effect model, and the results significantly reject the original hypothesis, which indicates that the fixed effect model based on within ordinary least square (OLS) estimation method is better than random effect model based on the feasible generalized least square (FGLS) estimation method. Furthermore, the modified Hausman statistic is used for a more accurate test, and the results also confirm the applicability of the fixed effect model. After testing, the panel had sequence correlation and heteroscedasticity, so the classical Hausman test was no longer applicable. Therefore, to address the heteroscedasticity problem, Hausman based on bootstrap method was used for further tests. The results further indicate that the fixed effect model should be chosen to obtain consistent estimation results.

Empirical Results
The results of the DID regression analysis are presented in Table 2. Columns (1), (3), (5), (7), (9) and (11) represent the results where we take the effect of meteorological factors on concentrations of air pollutants into consideration, that is, the meteorological factors are included in our empirical model as control variables, so as to obtain a more precise estimation of policy effect, while the other columns do not. Columns (1) and (2) show that the BSDW policy has significantly reduced the urban PM 2.5 concentration after its release. As shown in Column (1), the estimated coefficient of policy variable is −4.639 under 1% significance. This result is similar but larger than that of Jiang et al. [43], who found that PM 2.5 was reduced by 0.2637 µg/m 3 after the policy taking place. Column (2) represents that the estimated coefficient of policy variable is −1.951 under 1% significance. Obviously, there is no change in the direction and significance of the coefficient of policy variable. Wang et al. [73] showed that the main sources of PM 2.5 in BTH region were secondary nitrate (36-58%), traffic (8-26%), residential coal combustion (8-16%) and biomass burning (5-16%). Among them, secondary nitrate is formed by photochemical transformation of NO 2 with significantly high photochemical activity [74]. In line with the main sources of PM 2.5 , BSDW policy attaches great importance to implementing total coal consumption control and actively adjusting transportation structure and developing green transportation system. Therefore, the underlying reason for the emission reduction effect of BSDW policy lies in that it inhibits the production of pollutants from the source. Similarly, from the perspective of PM 10 (columns (3), (4)), SO 2 (columns (5), (6)), CO (columns (9), (10)) and O 3 (columns (11), (12)), the results also signify that the BSDW policy has significantly reduced the pollution in cities. As presented in column (10), the estimated coefficient of the policy variable is −0.057 under 1% significance. Column (9) controls the meteorological factors, and the results show that the corresponding coefficient is −0.014 under 1% significance. Similarly, as for O 3 , the estimated coefficient of the policy variable is −7.882 under 1% significance, which is shown in column (12). After controlling meteorological factors, the results are basically consistent. Notes: Robust standard errors are clustered at city level in parentheses; **, *** represent 5%, and 1% significance levels, respectively.
However, the concentration of NO 2 increased significantly when the BSDW policy came into force. As shown in column (7), the estimated coefficient of the policy variable after controlling the meteorological factors is 1.865 under 1% significance. Evidence exists that NO 2 has negative environmental-human health consequences, especially in circumstances where prolonged exposure is experienced [75]. Nitrogen oxides are mainly produced by anthropogenic activities, including vehicular emissions, biomass burning, fossil fuel combustion, thermal power plants, domestic solid fuel use, while natural sources include lightning and soils [76][77][78]. Matandirotya and Burger [79] found that seasonal biomass burning, which increased during the winter season, was also related to the emission of nitrogen oxides. However, the BSDW policy incentives the usage, promotion, and transformation of biomass power generation; specifically, it proposes to encourage the development of county biomass cogeneration, biomass briquette boilers and biological natural gas where resources are available. In view of the fact that the policy did not refine the specific indicators, the proposal to implement ultra-low emission transformation of biomass boilers may not reduce the emission of nitrogen oxides. Moreover, the efforts to reduce the emission of nitrogen oxides caused by domestic solid fuel use and thermal power plants may fail due to the lack of explicit emission reduction targets, although existing technologies such as selective catalytic reduction of NO x with NH 3 can make NO x emissions meet the standard [80].
The results of the DID regression analysis on AQI and different regions are presented in Tables A2 and A3. Based on the results, the BSDW has significantly reduced the urban AQI concentration and air pollutants in Beijing, Tianjin, and Hebei province, which shows that the results in this paper are robust.

Parallel Trend
We investigate whether the DID parallel trend is likely to exist. The original hypothesis of the parallel trend test is that before the implementation of the policy, the change trend of target variables over time is basically the same in the treatment group and the control group. To assess the plausibility of the parallel trend, the mean of the dependant variables for the treatment group and the control group can be plotted to assess the plausibility of the parallel trend [66,81]. Figure 2. indicates that the different groups appear to evolve along a similar temporal path, at least prior to treatment. After the implementation of the BSDW policy in the treatment group, a significant difference in slope can be detected after the implementation of the BSDW policy (in Figure 2 a,c,d,f). Thus, existing evidence indicates that parallel trends roughly exist, therefore, it can be concluded that the results based on DID method in Section 4.2 are robust.
variables for the treatment group and the control group can be plotted to assess the plausibility of the parallel trend [66,81]. Figure 2. indicates that the different groups appear to evolve along a similar temporal path, at least prior to treatment. After the implementation of the BSDW policy in the treatment group, a significant difference in slope can be detected after the implementation of the BSDW policy (in Figure 2 a,c,d,f). Thus, existing evidence indicates that parallel trends roughly exist, therefore, it can be concluded that the results based on DID method in Section 4.2 are robust.

Effect of Unobserved Variables
In addition to the parallel trend assumption test, we also test the effect of unobserved variables. The DID method is helpful for reducing the number of control factors to be considered in the model. To eliminate the effect of unobserved variables, this paper controlled the regional characteristics and adopted the two-way fixed effect model in the basic regression. However, it may still be difficult to observe and control some factors that change over time and space. Specifically, different industrial policies carried out by local governments in different regions may affect the upgrading of manufacturing industry, finally resulting in the estimation error. In order to eliminate the effect of unobserved factors, following the studies of La Ferrara et al. [82] and Liu and Lu [83], this paper uses the following method to indirectly test whether the unobserved characteristics of these regions will affect the estimation results.
According to the estimation in Equation (1), the formula for calculating the estimated coefficient of timei,t × treati,t is as follows:

Effect of Unobserved Variables
In addition to the parallel trend assumption test, we also test the effect of unobserved variables. The DID method is helpful for reducing the number of control factors to be considered in the model. To eliminate the effect of unobserved variables, this paper controlled the regional characteristics and adopted the two-way fixed effect model in the basic regression. However, it may still be difficult to observe and control some factors that change over time and space. Specifically, different industrial policies carried out by local governments in different regions may affect the upgrading of manufacturing industry, finally resulting in the estimation error. In order to eliminate the effect of unobserved factors, following the studies of La Ferrara et al. [82] and Liu and Lu [83], this paper uses the following method to indirectly test whether the unobserved characteristics of these regions will affect the estimation results.
According to the estimation in Equation (1), the formula for calculating the estimated coefficient of time i,t × treat i,t is as follows: In Equation (4), time i,t × treat i,t notes the interaction item of two dummy variables (treat i,t and time i,t ). β 3 (the coefficient of time i,t × treat i,t ) reflects the net effect of the policy. X i,t represents all the control variables included in this paper. i and t represent the monitoring station i and day t, respectively. ε i,t denotes random disturbance term. Theoretically, if γ = 0, it means that the unobserved factors will not interfere with the estimated results and thus β 3 is unbiased. However, it is difficult to directly test whether γ is 0. Therefore, if we replace the time i,t × treat i,t with a variable, which will not have real effect on the Air quality i,t in theory, and its coefficient is estimated to be zero, then γ = 0 can be derived inversely from this result. Thus, this paper randomizes the impact of the BSDW policy on specific regions (generated by computer), and then repeats this exercise 100 times. Such random processing can ensure that the BSDW policy will not affect the Air quality i,t , namely γ = 0. In this case, the mean value of γ is estimated, and the distribution of γ is shown in Figure 3. The average value of γ is close to zero and not significant compared with the basic results. It can be further found from Figure 3. that γ obtained from 100 random processes is distributed around 0, so we can deduce from the inverse that γ = 0. These results further boost our confidence that our findings are not severely biased by the unobserved time-variant regional characteristics.

Different Industry Type of Cities
Regional development difference is one of the prominent issues in China's sustainable and coordinated development [84]. Differences exist in the resource endowment and development path of different regions, leading to the industrial structures various in different regions [53]. Because the industrial structure has a decisive impact on resource consumption and pollutant emissions [85,86] and the content of the BSDW policy implemented in all cities is the same, we further classify cities according to industry types to estimate the emission reduction effect of the BSDW policy. Following the research of Yu et al. [44], this paper divides the treated cities in BTH region into two groups: one is urban consumption city (including Baoding, Beijing, Tianjin, Langfang, Qinhuangdao and Shijiazhuang); the other is resource intensive city (including Xingtai, Handan, Chengde, Zhangjiakou, Tangshan, Cangzhou and Hengshui). In accordance with the energy utilization characteristics of BTH region, including energy intensity and per capita consumption of implied energy, the cities in BTH region can be divided into urban consumption cities and resource intensive cities. Specifically, the urban consumption cities refer to the developed metropolis dominated by service industry and high-end manufacturing, with rapid economic development, low energy intensity but high implied energy per capita consumption. Resource intensive cities rely on local natural resources for development, and thus are energy intensive, with high energy intensity but low implied energy per capita consumption. They often rely on abundant energy to develop high-energy consumption sectors, such as mining sector. The proportion of industry is relatively high. The regression results are shown in Table 3.
times. Such random processing can ensure that the BSDW policy will not affect the Air qualityi,t, namely γ = 0. In this case, the mean value of γ is estimated, and the distribution of γ is shown in Figure 3. The average value of γ is close to zero and not significant compared with the basic results. It can be further found from Figure 3. that γ obtained from 100 random processes is distributed around 0, so we can deduce from the inverse that γ = 0. These results further boost our confidence that our findings are not severely biased by the unobserved time-variant regional characteristics.   Notes: Robust standard errors are clustered at city level in parentheses; **, *** represent 5%, and 1% significance levels, respectively.
In urban consumption cities, rows (1)- (6) show that the emission reduction effects of the BSDW policy on PM 2.5 , PM 10 , SO 2 , NO 2 , and O 3 are −1.751, −4.503, −1.634, 0.943 and −6.950 under 1% significance. However, the emission reduction effect of CO is no longer significant. In resource intensive cities, the results show that the emission reduction effects of the BSDW policy variable are −2.579, −3.487, −2.162, 3.021, −0.033 and −8.389, respectively. These results are consistent with the results in Section 4.2. However, it is obvious that the policy effect varies from region to region. In resource intensive cities, the BSDW policy can better play its emission reduction effect on PM 2.5 , SO 2 , CO and O 3 .
There are many energy-production and processing industries in resource intensive areas, such as mining and metallurgy [44]. Therefore, resource intensive cities have great potential for energy conservation, and should be the key area for energy conservation and emission reduction. Although there are abundant energy resources in the local area, the energy exploited has been encouraged to be used locally, processed into high-energy consumption products, and exported. Additionally, the energy cost and access threshold of high-energy consumption sectors are very low, and even encouraged. All of the above realities lead to the fact that resource intensive cities do not play the potential of energy conservation. On the one hand, the demand for high-energy consumption areas continues to grow. On the other hand, expanding the capacity of high-energy consumption sectors is the most convenient way for energy-oriented regions to promote economic growth, which drives the transfer of high-energy consumption sectors between the two types of regions, as well as the transfer of energy conservation and emission reduction pressure. In other words, the demand for high-energy consumption products in other regions is the reason for the pressure of energy conservation and emission reduction in resource intensive cities. The BSDW policy emphasizes that cities should speed up the adjustment of energy structure to build a clean, low-carbon and efficient energy system. Thus, resource intensive cities have great potential for energy transformation.

Key Environmental Protection Cities and Other Cities
Population and capital are most concentrated in cities, so is the environmental pressure. Speeding up the process of urbanization is an important way for China to achieve the goal of building a moderately prosperous society in an all-round way. China's urbanization is in the stage of accelerated development, urban environmental protection work will directly affect the success of national environmental protection work. On 26 November 2007, the State Council issued the "Eleventh Five Year Plan for national environmental protection" [87], which clearly stipulated 113 key urban environmental protection cities. The comprehensive prevention and control of air pollution in 113 key environmental protection cities such as Beijing and Tianjin are the focus of China's air pollution prevention, as well as the control and efforts to improve the quality of urban and regional air environment. The population of these 113 key environmental protection cities accounts for 50.1% of China's urban population, and the GDP accounts for 71.3% [43]. At the same time, the plan also clearly proposed to focus on 113 cities to comprehensively promote urban environmental protection work. Therefore, according to the plan, the treated cities in BTH region can be divided into two groups, one is the key environmental protection city, the other is not the key environmental protection city. The regression results are shown in Table 4. Rows (1)- (12) can further verify the robustness of the basic regression results. Whether it is a key environmental protection city or not, the BSDW policy will significantly reduce the concentration of air pollutants, except for NO 2 . For PM 10 , the estimated coefficient of the policy variable in key environmental protection cities is −4.323 under 1% significance, and the corresponding coefficient in other cities is −2.396 under 1% significance. As shown in Table 4, compared with the key environmental protection cities, the emission reduction effect of the BSDW policy in other cities is generally low. Figure 4 depicts the distribution of key environmental protection cities and other cities. The industrial distribution in BTH region plays an important role. Heavy industry leads to massive energy consumption (most of which is still provided by coal), and the most polluting is in the south of BTH region, which leads to serious air pollution [88].
Thus, the SO2 emission reduction effect in other cities is better than key environmental protection cities.

Heating Seasons
Existing literature [93][94][95] identified that heating in winter increased regional air pollution, which has become the focus of the BSDW policy [88]. Therefore, this paper further identifies the effect of the BSDW policy in heating seasons. According to the planning for heating seasons of the government, the heating season of the BTH region in this paper consists of two periods: the first period is from 1 January 2018 to 14 March 2018, and the second period is from 15 November 2018 to 31 December 2018. The results are displayed in Table 5. The results in columns (3), (5) and (6) are consistent with the basic regression results, that imply, after controlling the meteorological factors, the BSDW policy will significantly reduce the concentration of SO2, CO and O3 in heating seasons in the BTH region. For SO2, the estimated coefficient of the policy variable is −1.316 under 1% significance. As for CO According to the results in Table 3, the estimated coefficient of the BSDW policy in key environmental protection cities are always lower than those of other cities, except for SO 2 , CO and O 3 . This shows that after the implementation of the BSDW policy, compared with other cities, the emission reduction effect on PM 2.5 , PM 10 and NO 2 in key environmental protection cities is better. Therefore, we can consider including more cities in the list of key environmental protection cities to further improve the air quality in other cities. For SO 2 , however, some studies have pointed out that the combustion of coal in boilers is associated with the release of air pollutants, especially SO 2 and total suspended particles (TSP) [89]. The BSDW policy puts forward that "cities should further eliminate the combustion of coal in boilers". By comparing the government work reports of all cities in 2019, it can be found that almost all key environmental protection cities (such as Beijing, Tianjin, etc.) do not clearly mention the index of coal-fired boiler transformation, while most other cities (such as Xingtai, Zhangjiakou, etc.) have explicitly mentioned the results of coalfired boiler transformation. For example, Xingtai pointed out that 59 coal-fired boilers with 35 steam tons and below, 671 low nitrogen combustion boilers were eliminated [90]. Besides, 1663 coal-fired boilers were eliminated in Zhangjiakou [91] and "2695 coal-fired boilers were eliminated and banned" in Chengde [92] in the government word reports. Thus, the SO 2 emission reduction effect in other cities is better than key environmental protection cities.

Heating Seasons
Existing literature [93][94][95] identified that heating in winter increased regional air pollution, which has become the focus of the BSDW policy [88]. Therefore, this paper further identifies the effect of the BSDW policy in heating seasons. According to the planning for heating seasons of the government, the heating season of the BTH region in this paper consists of two periods: the first period is from 1 January 2018 to 14 March 2018, and the second period is from 15 November 2018 to 31 December 2018. The results are displayed in Table 5. Note: Standard errors are clustered by city level and reported in parentheses; ***, **, represent significant levels of 1%, 5%, respectively.
The results in columns (3), (5) and (6) are consistent with the basic regression results, that imply, after controlling the meteorological factors, the BSDW policy will significantly reduce the concentration of SO 2 , CO and O 3 in heating seasons in the BTH region. For SO 2 , the estimated coefficient of the policy variable is −1.316 under 1% significance. As for CO and O 3 , a significant decrease also exits, with the estimated coefficients of the policy variable being −0.023 and −8.558 under 1% significance, respectively. The most striking observation to emerge from the result is, as shown in columns (1), (2) and (4), this paper finds that the implementation of the BSDW policy does not reduce the concentration of PM 2.5 , PM 10 and NO 2 ; on the contrary, it causes a significant increase in the concentration of these air pollutants. For PM 2.5 , PM 10 and NO 2 , the estimated coefficients of the policy variable are 2.099, 15.126 and 5.527 under 1% significance, respectively. This is a rather disappointing result, which indicates that with the implementation of the BSDW policy, the concentrations of PM 2.5 , PM 10 and NO 2 increase significantly. However, this result has not previously been described. A possible explanation for these results may be since the lack of adequate natural gas supply sources, the original coal-fired heating facilities were not demolished, which induces that it is difficult to promote the "continue to promote coal to electricity and coal to gas" proposed by the BSDW policy.

City Heterogeneity
Since 2013, the emergency response plan of heavy air pollution in the BTH region has gradually tended to a unified emergency early warning response standard. Under different levels of heavy pollution weather warning, different cities and municipalities should formulate and implement cost-effective and differentiated heavy pollution emergency plans according to their actual situations, namely "one city, one policy" [96]. Therefore, this paper explores the different influences of the BSDW policy on different cities. The estimation results of the BSDW policy on different cities are shown in Table 6.
As shown in Figure 5, for PM 2.5 , the cities with the best pollutant reduction effect are Xingtai, Shijiazhuang and Handan; for PM 10 , Shijiazhuang, Beijing and Xingtai are the cities with the largest decrease in pollutant concentration; for SO 2 , emission reduction effect in Shijiazhuang, Xingtai and Handan are the best; for NO 2 , only Beijing has achieved significant emission reduction of NO 2 ; the concentration of CO decreases significantly in Hengshui, Cangzhou, Tangshan, Xingtai and Shijiazhuang; for O 3 , the pollutant emission reduction effect in Handan, Hengshui and Xingtai are the best. On the whole, the emission reduction effect of cities located in the south of the BTH region are better and the cities in the north of the BTH region have great potential for pollutant emission reduction. Based on the above results, the local governments of different cities should formulate specific measures to implement the BSDW policy according to their own characteristics.   Notes: Robust standard errors are clustered at city level in parentheses; *** represent 1% significance levels, respectively.
As shown in Figure 5, for PM2.5, the cities with the best pollutant reduction effect are Xingtai, Shijiazhuang and Handan; for PM10, Shijiazhuang, Beijing and Xingtai are the cities with the largest decrease in pollutant concentration; for SO2, emission reduction effect in Shijiazhuang, Xingtai and Handan are the best; for NO2, only Beijing has achieved significant emission reduction of NO2; the concentration of CO decreases significantly in Hengshui, Cangzhou, Tangshan, Xingtai and Shijiazhuang; for O3, the pollutant emission reduction effect in Handan, Hengshui and Xingtai are the best. On the whole, the emission reduction effect of cities located in the south of the BTH region are better and the cities in the north of the BTH region have great potential for pollutant emission reduction. Based on the above results, the local governments of different cities should formulate specific measures to implement the BSDW policy according to their own characteristics.

Conclusions and Policy Implications
There is an urgent need to improve air quality in China's BTH region. The three-year Action Plan to Win the BSDW has been implemented on a large scale. This paper has taken the implementation of the BSDW policy in the BTH region as a quasi-natural experiment and regarded the BSDW policy as an exogenous policy impact. Using daily panel data from 947 grid monitoring points in the BTH region in 2018, this paper sets out to analyze the emission reduction effect of China's air governance policies in BTH region by applying DID model. Furthermore, focusing on various pollutants, different city types, heating seasons and city level, this paper explores the different effects of the BSDW policy. The main conclusions are as follows.
Our findings propose that the BSDW policy has partly significantly improved overall air quality and reduced atmospheric pollutant emissions in the BTH region. The policy leads to a reduction in the daily concentration of PM 2.5 , PM 10  The parallel trend test results and the test of effect of unobserved variables indicate that the reduction in air pollution concentration after 4 July 2018 is due to the BSDW policy, rather than other factors. Substituting the dependent variable with AQI and dividing the treatment group into three groups obtained consistent results; the concentration of AQI (−3.012) has decreased significantly after controlling the meteorological factors, indicating that the conclusions in this paper are robust. The heterogeneity of different city types and different cities also shows that the aforementioned conclusions are robust. On the whole, the emission reduction effect on cities located in the south of the BTH region is better. Resource intensive cities, non-key environmental protection cities and cities in the north of the BTH region have great potential for pollutant emission reduction. The results of the heterogeneity analysis show that the estimated coefficients of the policy variable in key environmental protection cities are always lower than those of other cities, except for SO 2 , CO and O 3 .
Despite these promising results, questions remain. One of the more significant findings to emerge from this study is that with the implementation of the BSDW policy, the concentration of NO 2 increases significantly, with the estimated coefficient of 0.673 under 1% significance. The corresponding coefficient is 1.865 under 1% significance after controlling the meteorological factors. Another important finding is that in heating seasons, the BSDW policy will reduce the concentration of SO 2 (−1.316), CO (−0.023) and O 3 (−8.558) under 1% significance. Surprisingly, the concentrations of PM 2.5 (2.099), PM 10 (15.126) and NO 2 (5.527) increase significantly in heating seasons, which might be of interest to policymaking.
According to the findings, several suggestions have been proposed to further promote the implementation of the BSDW policy. First, set and aim to achieve a specific goal based on city level to expedite an optimal "Blue Sky" strategy. Different cities should formulate and implement differentiated air pollution prevention plans according to local actual situations, namely "one city, one policy". In addition, incorporating more cities into the list of key environmental protection cities may contribute to the improvement of the air quality in other cities. Second, refine the specific emission reduction indicators related to domestic solid fuel use and thermal power plants. Besides, attention should be paid to accelerate the transformation of gas-fired boilers with low nitrogen and urban biomass boilers with ultra-low emissions and formulate specific transformation indicators. Third, further improve the natural gas supply system in the BTH region, and clean heating diffusion should be encouraged and enforced. Data Availability Statement: The data are not publicly available due to privacy.

Conflicts of Interest:
The authors declare no conflict of interest. Appendix A Figure A1. Spatial distribution of sample monitoring points. 100 160 160 120 --100 Figure A1. Spatial distribution of sample monitoring points.  Notes: Robust standard errors are clustered at city level in parentheses; **, *** represent 5%, and 1% significance levels, respectively.

Appendix B
Step 1.