This study refers to Chen et al. (2013) and uses a quasi-natural experiment on centralized heating policies in winter in northern China (north of the Qinling-Huaihe River boundary) to address the possible endogeneity problem between air pollution and firm productivity using regression discontinuity and Stata software. The reason for choosing the Qinling-Huaihe River is that there is a significant difference in the average January temperature between the south and north areas of the boundary; therefore, the centralized heating policy in winter in the north is an exogenous policy in a quasi-natural experiment. This study sets the following estimating equation. D i t = 1 ,   L c ≥ 0 0 ,   L c < 0 (1) T F P i t = δ 0 + δ 1 D i c t + δ 2 f L c + ∑ n = 2 k δ n * C o n i t + τ i + γ t + ε i t (2) P M 2.5 i t = γ 0 + γ 1 D i c t + γ 2 f L c + ∑ n = 2 k γ n * C o n i t + τ i + γ t + ε i t (3) T F P i t = α 0 + α 1 P M 2.5 i t + α 2 f L c + ∑ n = 2 k α n * C o n i t + τ i + γ t + ε i t (4) where c denotes the city; the explanatory variable T F P i t represents the productivity of firm i in year t; P M 2.5 i t is the air quality of the city where the firm is located; C o n i t represents a set of control variables; τ i , γ t , and ε i t are the individual and time dummy variables and random disturbance terms, respectively; L c is the driving variable, which denotes the difference between the latitude of the city and the Qinling-Huai River boundary ² D denotes the indicator variable, which takes the Qinling-Huai River boundary under the background of China’s centralized heating policy as a discontinuity. When city c is located north of the boundary, that is, the city is located north of the Qinling-Huai River, L c >0 and D takes 1. Otherwise, when city c is located south of the boundary, that is, the city is located south of the Qinling-Huai River, L c <0 and D takes 0. f L c is the higher-order polynomial adjustment function of the driving variable L c .

Referring to Tang et al. (2017), there is a spatial autocorrelation effect of environmental pollution in cities. Therefore, considering that firm productivity is affected by local urban air pollution and also by the spatial spillover effects of air pollution from other cities, we use local pollution and pollution from surrounding areas as the core variables to explore the total impact of local firm productivity by local pollution and pollution spillover from surrounding areas. Through the Moran index, it can be shown that there is a spatial autocorrelation effect of air pollution PM2.5 in the city where the enterprise is located, and this is consistent with the conclusions obtained by Fan et al. (2019) through the spatial correlation of haze pollution measurement ³ Next, we add a spatial econometric model, as shown in Equations 5, 6. T F P i t = β W P M i t + β 0 + β 1 P M 2.5 i t + ∑ n = 2 k β n * C o n i t + τ i + γ t + ε i t (5) W P M i t = ∑ j = 1 c = 214 W i j * p m j t (6) where W P M i t represents the total impact of air pollution in the city of other enterprise j on the city of local enterprise i, that is, ∑ j = 1 c = 214 W i j * p m j t . W is the 214*214 spatial weight matrix, and this study uses MATLAB software to calculate the distance between cities and form the inverse geographic distance weight matrix W i j by using the latitude and longitude data of 214 major cities. W i j denotes the inverse geographic distance weight matrix of the city where firm i is located and the city where other firm j is located; p m j t denotes the air quality of the city where other firm j is located in year t, and the rest of the variables are consistent with the previous section.

The subject of this study is firm productivity, traditionally measured using ordinary least squares and fixed effects; however, these two methods may cause serious bias and endogeneity problems. There are two semi-parametric methods of OP and LP proposed by Olley and Pakes (1996) and Levinsohn and Petrin (2003) to estimate firm productivity. In the OP method, it is assumed that a manufacturer’s investment is a strictly increasing function of its own productivity, so the corresponding productivity can be determined by deriving the inverse function of the investment demand function. The OP method introduces survival probability and investment amount as proxy variables to solve the problems of selectivity and simultaneity bias, but the firm’s real investment value must be greater than zero. Thus, samples with zero investment are not estimated. However, as not every firm invests every year, many firm samples are discarded in the estimation process. Meanwhile, the LP method has lower data requirements. It does not use the amount of investment as a proxy variable and can instead use the input indicators of intermediate goods, which allows the researcher the flexibility to choose data. Therefore, we measure firm productivity using the LP method to determine the response variables in the benchmark regression.

The core explanatory variable in this study is air quality, as measured by the PM2.5. Air pollutants currently emitted in China are mainly sulfur dioxide, carbon dioxide, soot, and PM. Since the PM2.5 can synthetically include the above pollutants emitted from the combustion of fuels and further chemical reactions in the air, it is useful for the analysis of gaseous pollutants emitted from production enterprises and suitable for measuring air quality. Owing to the short time span and difficulty in obtaining pollution data in China, there may be problems such as data manipulation and falsification. Therefore, we used the world PM2.5 density map from 1998 to 2016, publicly available at the Data Center of the Astronaut Earth Observing System Data and Information System hosted by Columbia University ⁴ Then, we positioned each raster to municipal administrative units by ArcGIS software, and then averaged the raster data falling in the municipal administrative units as the value of PM2.5 emission concentration of cities. Unlike ground-observed air pollution data, satellite-extracted data are better able to match mobile population data in time and space while avoiding problems such as data manipulation and human measurement bias (Ghanem and Zhang 2014).

Based on the mechanism analysis, we select the number of patent applications, number of enterprise employees, and government subsidies as mediating variables to represent the enterprise R&D innovation effect, human capital effect, and policy spillover effect, respectively, to test the mechanism of the effect of air quality on enterprise productivity. Among them, enterprise patents include three types of invention patents, utility model patents, and design patents. We choose the sum of three types of patent applications plus one to take the logarithm to measure the number of enterprise patent applications (Hall and Harhoff 2012); to narrow the data differences and obtain smoother data, the number of enterprise employees and government subsidies are logarithmically processed in this study.

In addition to the required variables to estimate air pollution and firm productivity, we include the following control variables. (1) Firm size. As the size of the firm increases, the firm will have a higher scale effect, stronger finances, a more standardized management system, and usually higher firm productivity (Lucas 1978; Majumdar 1997). (2) Return on total assets. We choose the ratio of the firm’s net profit for the year and total assets at the end of the period to measure the return on total assets. The more profitable a firm is, the more likely it is that sufficient funds will be available for R&D innovation and high-quality production activities (Bogliacino and Pianta 2013). (3) Fixed assets ratio and capital expenditure ratio. Fixed assets are a category of assets formed when capital expenditure is capitalized. When the fixed assets ratio or capital expenditure ratio is lower, the liquidity of the enterprise’s assets is faster and the capital operating capacity is stronger (Wang 2002). (4) Current gearing ratio and gearing ratio. Both reflect a firm’s financing constraint ability to a large extent. When the current gearing ratio is lower, the firm may face higher financial risk, which reduces productivity (Ayyagari et al., 2010). (5) Cash flow. In general, cash flow has a positive effect on firm productivity, and firms maintain a certain amount of cash flow to avoid the shortage of funds (Chen and Guariglia 2013). (6) Firm age: There may be two effects of firm age on firm productivity. On the one hand, as firms age, they accumulate more technical information, beneficial to firm productivity (Farinas and Moreno 2000). On the other hand, new market entrants may be more dynamic; in this case, firm productivity decreases as firms age (Brandt et al., 2012). Therefore, we add firm age and the squared terms of firm age as control variables. The specific expressions are shown in Table 1.

Considering that the data on listed companies in China before 2007 are missing and/or inadequate, the sample period selected for this study is 2008–2016. Our sample includes A-share companies listed on the Shanghai and Shenzhen stock markets by securities code and year identifier. We processed the data by removing the following: (1) missing control variables, (2) enterprises established after the listing time, and (3) and financial and insurance services and enterprises containing ST, ST*, and PT categories ⁵ Then, the filtered data were matched with the PM2.5 level for 214 major cities. Next, most of the variables were logarithmically processed, resulting in a final sample of 8,594 observations from the 1,829 enterprises listed. Table 2 reports the descriptive statistics of all variables. To determine the viability of the variables, we performed multicollinearity and heteroscedasticity tests; the results showed that we did not need to worry about multicollinearity and homoscedasticity among the variables ⁶

Before conducting regression discontinuity estimation, it is first observed whether air pollution and firm productivity show discontinuity change and the results are shown in Figure 1. Figure 1 shows that regardless of whether the fit is linear or quadratic, both the productivity and PM2.5 of enterprises jump significantly at the discontinuity, that is, the PM2.5 emissions in the northern region are significantly higher than those in the southern region, this conclusion agrees with the findings of Chen et al. (2013). While the enterprises in the southern region have significantly higher productivity, indicating that relatively good air quality is beneficial to enterprise productivity. In addition, to accurately identify the causal relationship between variables, we also perform an endogeneity test in the 5.3 section.

The results of this estimation are reported in Table 3; column (1) shows that the central heating policy has a significant impact on enterprise productivity, that is, the central heating policy makes the enterprise productivity in the southern region significantly higher than that in the northern region at the level of 1%. When control variables and fixed effects are added in columns (2) and (3), the regression results do not change significantly. Similarly, the results in column (4) show that the central heating policy has a significant impact on air pollution. After introducing control variables and controlling for fixed effects in columns (5) and (6), the central heating policy makes the average PM2.5 concentration in the southern region about 0.23 µg/m³ lower than that in the northern region. The estimation results in column (7) show that air pollution significantly reduces the productivity of enterprises. Columns (8) and (9) gradually add the corresponding control variables and fixed effects, and the regression results are all significant at the 1% level, that is, with every 1% increase in PM2.5 concentration, firm productivity decreases by 0.32%.

In addition, we consider the air quality index (AQI) as a surrogate indicator of air pollution PM2.5 for robustness testing, which monitors six pollutants: sulfur dioxide, nitrogen dioxide, PM10, PM2.5, carbon monoxide, and ozone. We take the natural logarithm of AQI, and the regression results are reported in column (10) of Table 3. The results show that after adding a series of control variables and fixed effects, AQI still significantly reduces enterprise productivity, which is consistent with the explanatory variable PM2.5, showing the robustness and reliability of the conclusions.

Considering that firm productivity may be affected by its own lags, the model is constructed to examine its dynamic effects, as expressed in Equation 7. T F P i t = β 0 + β 1 T F P i , t − 1 + β 2 P M 2.5 i t + ∑ n = 3 k β n * C o n i t + τ i + γ t + ε i t . (7)

Eq. 7 indicates that the productivity of firm i in year t-1 will affect that in year t. As indicated in column (1) in Table 4, the coefficient of PM2.5 is −0.01, which is highly significant. In addition, the coefficient of TFP_i,t-1 is significantly positive. This result implies that the higher a company’s productivity in the previous year, the higher its productivity in the second year; that is, it has a certain “inertia.”

Considering unobservable individual effects, this may lead to inconsistent estimation results and endogeneity problems. In dynamic regression, GMM tests are generally used to solve the endogeneity problem, but compared to the systematic GMM, the instrumental variables of differential GMM are often weak. Therefore, we use the systematic GMM method to re-estimate this model to solve the endogenous problem of fixed effects (Blundell and Bond 1998), as presented in column (2) of Table 4. The regression results of the systematic GMM show that the first-order lagged terms of the explanatory variables pass the significance test, and this result is consistent with the fixed effects regression results, which indicate that an increase in PM2.5 effectively reduces firm productivity, further strengthening the reliability of our results.

In addition, at the theoretical level, the systematic GMM test needs to satisfy two hypotheses. First, the validity of the instrumental variables used is tested by the Sargan or Hansen overidentification constraint test, based on the original assumption that the instrumental variables are uncorrelated with the error term. The p-value of the Hansen test is greater than 0.1 in column (2) of Table 4, indicating that the null hypothesis is accepted; that is, the instrumental variables of the systematic GMM model are valid. Second, the second-order serial correlation of the random error terms of the difference equation is tested by the Arellano-Bond autocorrelation test, where the original hypothesis is that there is no second-order serial correlation in the random error terms of the first-order difference equation. If the original hypothesis is not rejected, it implies that the instrumental variables are valid and the model is set up correctly. Therefore, we need to focus on the p-value of the second order autoregressive model (AR (2)). The p-value of the AR (2) test in column (2) of Table 4 is greater than 0.1, indicating that the null hypothesis is accepted; that is, there is no second-order autocorrelation problem in the residual term of the systematic GMM model difference, and the model setting is reasonable.

Owing to the global financial crisis in 2008, China’s economic growth rate has dropped significantly since 2013 and has entered a new normal mode. Therefore, we reselect the data from 2008 to 2012 and from 2013 to 2016 as samples to eliminate exogenous effects and shocks caused by the financial crisis and new economic normal. Table 4 presents the regression results, with air pollution still having a significant inhibitory effect on corporate productivity, consistent with previous empirical findings.

Table 5 shows the selection process of the spatial effects model and the corresponding regression results. First, from the log-likelihood function values, the SDM has greater values than the spatial lagged model (SLM) and the spatial error model (SEM), which indicates that the SDM model is the best fit. Second, we test whether the SDM model degrades to the SLM or SEM model, and the likelihood ratio (LR) and Wald statistics both reject the original hypothesis at the 1% level, indicating that the SDM model cannot degenerate into an SLM or SEM model. Therefore, the SDM is determined to be most suitable for this study.

Based on the existence of spatial autocorrelation of pollution in the city where the firm is located, column (3) of Table 5 reports the regression results of PM2.5 in the city where the firm is located and pollution spillover from surrounding cities on firm productivity under the SDM. The results show that both PM2.5 in the firm’s city and pollution spillover from surrounding cities significantly reduce firm productivity, indicating that firm productivity is not only negatively affected by local air pollution but also by PM2.5 spillover from surrounding cities.

China’s regional development is unbalanced and the terrain system is extremely complex; we should not ignore the differences in different regions. Therefore, this paper divides enterprises into 11 eastern enterprises such as Beijing, eight central enterprises such as Shanxi Province, and 12 western enterprises such as Sichuan Province according to the provinces where the enterprises are located ⁷ and investigates the impact of PM2.5 on enterprise productivity in different regions. Meanwhile, purely distance-based spatial weights cannot effectively capture spatial relationships, especially if the phenomenon is environmental in nature. Therefore, in addition to the inverse distance spatial weight matrix W1 calculated by latitude and longitude, we also select the economic distance type spatial weight matrix W2 constructed by urban per capita GDP and the economic-geographic distance type spatial weight matrix W3, constructed by multiplying the two and normalizing the result.

According to the regression results in Table 6, under the three spatial weight matrix settings, the productivity of enterprises in the eastern, central, and western regions will be negatively affected by PM2.5. Poor air pollution will have a greater negative effect on the productivity of eastern enterprises, indicating that enterprises in eastern provinces are more sensitive to air pollution. In addition, the WPM2.5 coefficient is significantly negative at the levels of 1%, 5%, and 10%, indicating a significant spatial correlation between air pollution at similar geographical locations or similar levels of economic development. Pollution spillover from surrounding areas has a negative spillover effect on the productivity of enterprises in the central area, and enterprises in the central and western regions are more affected by PM2.5 in surrounding cities.

Although regression discontinuity can effectively solve the endogeneity problem, it can only be locally random (near the discontinuity), resulting in weak external validity. We have to further consider the endogeneity of the air pollution question. Specifically, on the one hand, air pollution may drag down enterprise productivity by affecting the human capital, government subsidies, and R&D innovation capabilities; on the other hand, enterprise productivity itself may also affect air pollution levels through factors such as economic activities and population aggregation. Therefore, in order to solve the endogeneity problem, referring to Broner et al. (2012); Hering and Poncet (2014), we use the air flow coefficient (Airflow) as an instrumental variable of air pollution. The reason for this is that air mobility is correlated with air pollution, and the air mobility coefficient depends on wind speed and atmospheric boundary layer height, which is not directly related to enterprise productivity, thus satisfying the validity and exogenous assumptions of instrumental variables (Broner et al., 2012; Hering and Poncet 2014). Meanwhile, we use GS2SLS to effectively control for endogeneity and spatial spillover effects.

The regression results are shown in Table 7. We still conduct inspections from the east, middle, and west, and report the GS2SLS two-stage estimation results based on the economic-geographical weight matrix (W3) ⁸ . All the regressions control for the control variables presented above and the fixed effects of year and city. According to the first-stage regression results, air flow is significantly negatively correlated with PM2.5 concentration, indicating that stronger air flow is more conducive to alleviating air pollution. The LM statistic and p-value at the 1% significance level rejected the null hypothesis of “no endogeneity in explanatory variables,” indicating that air pollution is an endogenous variable and that the endogeneity issue must be addressed in this paper. Furthermore, the F statistic is much greater than 10, indicating that air flow and PM2.5 are highly correlated, and thus rejecting the null hypothesis of “instrumental variables are weak instrumental variables.”

In the regression results of the second stage, PM2.5 significantly inhibited enterprise productivity. In terms of spatial latitude, the coefficient estimation result of the spatial lag variable WPM2.5 based on the nested weight matrix W3 of economic-economic-geographical distance is significantly negative at the levels of 1%, 5%, and 10%, indicating that air pollution in surrounding cities has a significant negative spatial spillover effect on the productivity of firms located in central cities. It is consistent with the basic regression results, indicating that the conclusions are still reliable after dealing with the endogeneity issue.

Different air pollutants may have different impacts on firm productivity. In this study, we match existing data with the newly released environmental research data of listed companies in the CSMAR database to obtain information such as emissions of listed companies by subdivision. However, after distinguishing the subsamples of two air pollutants, sulfur dioxide and soot emissions, only a small amount of data remains, as shown in columns (1) and (2) of Table 8. The regression results show that the coefficients of both air pollutants are negative, but neither is significant. To compensate for the above lack of sample size, this study collects sulfur dioxide and soot emissions from each city to match the existing data, as shown in columns (3) and (4) of Table 8. From the regression results, it can be seen that the data collected from each city and the data from the CSMAR database of environmental studies of listed companies are consistent, that is, the coefficients of sulfur dioxide and soot emissions are both negative, but neither is significant. This indicates that these two air pollutants bring much less inhibitory effect on firm productivity compared to PM2.5, which reflects the importance of choosing PM2.5 as the core variable in this study.

By the end of 2020, the total assets and operating income of state-owned enterprises (SOEs) in China reached 218.3 and 59.5 trillion yuan, respectively, and their average annual growth rate during the “13th Five-Year Plan” period was 12.7% and 7.4%, respectively. After years of reform and development, SOEs have laid a solid foundation for China’s new stage of development. Simultaneously, the rapid development of the non-public sector has played a pivotal role in labor employment, technological innovation, national taxation, and economic development. Compared with SOEs, the vitality of non-SOEs is more affected by external policies and the environment. Meanwhile, compared with non-SOEs, the factor market price distortion of SOEs is more serious. Therefore, what is the difference in the impact of air quality on the productivity of SOEs and non-SOEs? For a better comparative analysis, columns (1) and (2) of Table 9 report the results of these two types of samples. The results show that the coefficient of PM2.5 is significantly negative at the 1% level, indicating that air pollution reduces productivity in SOEs and non-SOEs, with a greater impact on SOEs. Possibly, in China, the central and local governments control or invest in SOEs, which are not as sensitive as non-SOEs with respect to the cost pressure caused by air pollution control and the green innovations required to improve air quality. In addition, compared with non-SOEs, SOEs have insufficient information on the improvement of air quality, since the free flow of production factors is not flexible enough and the efficiency of resource reallocation is not high enough, which causes air pollution to have a greater impact on the productivity of SOEs.

Further, innovation is the first driving force for development. Generally, high-tech enterprises have a strong sense of innovation and carry out R&D activities continuously to enhance their core independent intellectual property rights. Columns (3) and (4) of Table 9 report the regression results of dividing the total sample into high-tech and non-high-tech companies, according to their industry category. The results indicate that the coefficients of PM2.5 are all significantly negative. Specifically, the former coefficient is −0.014, and the absolute value is smaller than the latter’s coefficient of −0.031. This shows that air pollution reduces the productivity of enterprises in all industries. Compared with enterprises in high-tech industries, air pollution has a higher negative impact on enterprises in non-high-tech industries, which is contrary to the research of Cao et al. (2022). One possible reason is that human capital capabilities and contribution to corporate productivity of high-tech companies are greater. Therefore, in the face of air pollution, companies will try to make up for the environmental requirements of high-quality labor, and the productivity of high-tech enterprises will not be hit as badly as that of non-high-tech enterprises.