According to Zhu et al. (2018), time series econometric techniques provide a feasible research method for identifying the asymmetric interaction effects of haze pollution. Under the time series econometric analysis framework, the asymmetric interaction effects of haze pollution can be identified using at least three tools, which are the traditional linear Granger causality test, nonlinear Granger causality test, and variance decomposition. Therefore, the research framework of asymmetric interaction effects of haze pollution and the construction of spatial weight is shown in Figure 1.

Linear Granger causality test is used to identify asymmetric spatial interaction. The change of haze pollution in a region may affect the change of pollution in other regions. Therefore, this paper identifies the asymmetric interaction effects of haze pollution under the vector autoregressive regression model (VAR) framework. In the VAR framework, the linear Granger causality test approach treats all variables as endogenous and does not require too many a priori constraints on the relationship between variables. It can reveal the asymmetric probability dependence between haze pollution time series and can be used to determine whether there is a lead-lag relationship between two series of haze pollution (Granger, 1969). If the previous information of one time series cannot help to improve the prediction ability of the current value of another time series, the former is said to be not the Granger cause of the latter; otherwise, the former is said to be the Granger cause of the latter. For example, suppose the current or previous information of air pollution in city X helps to improve the predictive ability of future values of haze pollution in city Y. In that case, X is the Granger cause of Y. This paper performs the Granger causality test according to the method proposed by Hsiao (1981). And with the help of the AIC criterion, we select the optimal lag order (Kilian, 2001).

The Nonlinear Granger causality test is used to identify asymmetric spatial interaction. As Granger and Newbold (2014) pointed out, “the real world is full of nonlinear relationships, and nonlinear models can help us simulate the real world more accurately.” For this reason, Diks and Panchenko (2006) proposed a nonlinear Granger causality test based on Hiemstra and Jones (1994). We draw on the idea of Diks and Panchenko (2006) to construct a bivariate VAR model based on a smooth time series. And we filter the linear components among the variables by extracting the VAR residuals and then conduct a nonlinear test on this basis. This paper uses the BDS method (Broock et al., 1996) to identify the air pollution spatial interaction by using the nonlinear Granger causality test under the premise of confirming the existence of a significant nonlinear dynamic trend between the series. In addition, we argue that if the original hypothesis of no nonlinear Granger causality is accepted for all lag orders (1–8), then there must be no asymmetric spatial interaction between haze pollution variables. If the original hypothesis of no nonlinear Granger causality is significantly rejected at all lags, we cannot rule out the existence of asymmetric spatial interaction among the haze pollution variables.

Variance decomposition is used to identify asymmetric spatial interaction. There are still some limitations in the nonlinear Granger causality test. Firstly, the nonlinear Granger causality test only considers the direction of the interaction effect of haze pollution but not the strength of the effect. Secondly, linear Granger and nonlinear Granger are based on a pure information flow perspective, which ignores the characteristics of haze pollution decaying with geographical distance. Whereas the variance decomposition technique, by analyzing the magnitude of the impact of each city’s haze pollution on the degree of another city’s haze pollution, provides relatively important information of all the stochastic perturbation that impacts the variables in the VAR model. Thus, this paper constructs a set of cities that may have asymmetric spatial interaction with each city, based on which VAR modeling is used to reveal the probability dependence and the degree of interaction effects of haze pollution among these cities. However, it is worth noting that the derivation of variance decomposition involves Cholesky decomposition in a certain order of variables, which makes the final analysis depends on the ranking of variables (Campbell, 1991) and cannot take into account the geographical distance factor. Therefore, in this paper, we take the idea of the first law of geography by Tobler (1970) into account. Due to the geographic decay law of the spatial interaction of haze pollution, the intensity of the interaction of haze pollution decreases with the increase of geographic distance. According to the above law, we get the basis for the ranking of variables, i.e., the closer the city is to a city in the Cholesky decomposition process, the more distant the city is, the more distant it is, the more distant it is. This way solves the problem of variable selection and variable ranking of variance decomposition in the VAR model and introduces the law of haze pollution decaying with geographic distance into the variance decomposition.

Anselin and Florax (2012) argued that the better the spatial weight matrix match the spatial structure of data, the better the fit of the model, and the better the explanatory power. In contrast, the existing weights, both symmetric and asymmetric weights, are given a priori and do not change over time, thus neglecting the endogeneity of the weight matrix (Harris et al., 2011; Corrado and Fingleton, 2012). Therefore, to more accurately identify the asymmetric spatial interaction and then construct spatially asymmetric weights that highly match the true spatial structure of real data. This paper identifies the asymmetric spatial interaction of haze pollution in the time series econometric analysis framework and constructs asymmetric spatial weights based on the identification results.

Spatial weight construction based on Granger causality test. In this paper, the spatial weights W₁ and W₂ are constructed according to the identification results of the linear Granger causality test and the nonlinear Granger causality test respectively, and the matrix elements are shown in Eq. 1. w i j = { 1 ( j → i ) 0 ( i = j , j → i ) (1) In Eq. 1, w _ij is a matrix element, i, j is the sample city (same below). j → i shows that the haze pollution in the jth city is the linear (nonlinear) Granger cause of the ith city, then w _ij = 1. j → i shows that the haze pollution in the jth city is not the linear (nonlinear) Granger cause of the ith city, then w _ij = 0; when i = j , w _ij = 0.

Spatial weight construction based on variance decomposition. The spatial weights W₃ are constructed according to the results of variance decomposition identification, and the matrix elements are shown in Eq. 2. w i j = [ 0 δ 12 ⋯ δ 1 n δ 21 0 ⋯ δ 2 n ⋮ ⋮ ⋱ ⋮ δ n 1 δ n 2 ⋯ 0 ] (2) Where, δ _ij is the contribution of haze pollution derived from variance decomposition, which indicates the degree of influence of city j on the city i. Further, the general form of its weight can be transformed into Eq. 3. w i j = { δ i j 0 ( j → i ) ( i = j , j → i ) (3) If at j → i , i.e., city j haze pollution has an impact on the city i haze pollution and the degree of influence is δ _ij , then w _ij = δ _ij ; if at j → i , i.e., city j haze pollution has no influence on the city i haze pollution, then w _ij = 0; and when i = j , w _ij = 0.

Given that the common theoretical framework for studying environmental pollution impact factors is the STIRPAT model proposed by York et al. (2003). Therefore, we use the STIRPAT model as the basis for selecting impact factors and the model is shown in Eq. 4. I i t = α P i t b A i t c T i t d μ i t (4) In Eq. 4, i and t denote city and year respectively, I _it represents environmental impact, P _it represents the population, A _it represents wealth, T _it represents technology, ε _it is the random error term. The STIRPAT is a multivariate nonlinear model, and the linear model is obtained by taking the logarithm of the equation at the same time. The STIRPAT model allows the introduction of other influential factors to analyze. According to the above analysis framework, this paper constructs a spatial econometric analysis model to investigate the social and economic impact factors of haze pollution by combining the spatially non-homogeneous characteristics of haze pollution. According to the existing research, the spatial econometric model should start with the spatial Durbin model. Because we can include the spatial lagged term of the explained variable (WY) and the spatial lagged term of the explanatory variable (WX) into the spatial Durbin model at the same time, which if ignored would result in a loss of consistency in the estimation of the remaining estimates (Elhorst, 2010; Elhorst, 2014). Given this, this paper develops an econometric model as in Eq. 5. Y t = ρ W Y t + α X t + β W X t + μ t (5) In Eq. 5, t is the year. The explanatory variable Y _t is haze pollution, X _t is the relevant influential factor and μ _t is the random error term. W _ij represents the asymmetric spatial weight matrix. ρWY _t and βWX _t represent the spatial dependence of the explanatory and explanatory variables respectively. ρ is the spatial lag coefficient, which reflects the influence of haze pollution in other cities on the sample city in the current period. β is also the spatial lag coefficient, which represents the influence of other urban impact factors on haze pollution in the sample city in the current period.

Studies have been conducted without spatial correlation tests when adding the spatial lagged terms of explanatory variables (WX) in the construction of spatial Durbin models (Hao and Liu, 2015; Wang and Fang, 2016; Du et al., 2018). Thus, adding all the spatial lagged terms of all explanatory variables to the model can lead to model over-explanation and problems such as multicollinearity. Given this, when adding the WX term in this paper, the Bivariate Moran test (Anselin, 2010) is performed first. Meanwhile, the regression coefficients of the above model cannot directly express the effect of explanatory variables on the explained variable, so further effect decomposition is needed. And we first rewrite the model as in Eq. 6. Y t = ( I − ρ W ) − 1 + ( I − ρ W ) − 1 ( α X t + β W X t ) + ( I − ρ W ) − 1 μ t (6) The matrix of partial derivatives of E (Y _t ) (the expected value of Y _t ) corresponding to the kth explanatory variable in X _t from city one to city n can be written as Eq. 9. These partial derivatives represent the effect of a one-unit change in a particular explanatory variable in one city on the explanatory variables in all other cities, which can be divided into direct effects and indirect effects (Mao et al., 2022). In Eq. 9, I _N is the unit matrix, d ¯ denotes the operator that calculates the mean of the diagonal elements of the matrix, and r s u m ¯ denotes the operator that calculates the row and means of the non-diagonal elements of the matrix (Elhorst, 2014). Direct effect and indirect effect are shown in Eqs 7, 8. Direct   effect:   [ ( I − ρ W ) − 1 ( α k I N + β k W ) ] d ¯ (7) Indirect   effect:   [ ( I − ρ W ) − 1 ( α k I N + β k W ) ] r s u m ¯ (8) Total   effect:   [ ( I − ρ W ) − 1 ( α k I N + β k W ) ] d ¯ + [ ( I − ρ W ) − 1 ( α k I N + β k W ) ] r s u m ¯ (9)

Global Spatial Autocorrelation is commonly used to analyze the distribution characteristics, which is generally portrayed by Moran’s I (Moran, 1950) and Geary’s C (Geary, 1954). Firstly, the Global Spatial Autocorrelation of urban haze pollution in China is examined using Moran’s I index, and the specific measurement is carried out according to Eq. 10. Where, S 2 = 1 / n ∑ i = 1 n ( x i − x ¯ ) 2 , x ¯ = 1 / n ∑ i = 1 n x i _, w _ij is the spatial weight matrix, n is the number of cities, and x _i denotes the observed value of haze pollution in the ith city. Mora n ′ s I = n ∑ i = 1 n ∑ j = 1 n w i j ( x i − x ¯ ) ( x j − x ¯ ) ∑ i = 1 n ∑ j = 1 n w i j ∑ i = 1 n ( x i − x ¯ ) 2 = ∑ i = 1 n ∑ j = 1 n w i j ( x i − x ¯ ) ( x j − x ¯ ) S 2 ∑ i = 1 n ∑ j = 1 n w i j (10) Next, Geary’s C index is used to test the Global Spatial Autocorrelation of urban haze pollution in China, measured according to Eq. 11. Unlike Moran’s I index, the value of Geary’s C index greater than one indicates a negative correlation, equal to one indicates no correlation, and less than one indicates a positive correlation. Geary ' C = ( n − 1 ) ∑ i = 1 n ∑ j = 1 n w i j ( x i − x j ) 2 2 ( ∑ i = 1 n ∑ j = 1 n w i j ) [ ∑ i = 1 n ( x i − x ¯ ) 2 ] (11) Using the above method, this paper measures the global Moran’s I and Geary’s C values of PM2.5 in 160 cities. Moran’s I is significantly positive under all three weighting schemes, while Geary’s C is significantly less than one for all of them. This indicates a significant positive spatial correlation characteristic of haze pollution in 160 cities. Therefore, the spatial characteristics of haze pollution should be fully considered; otherwise, it will lead to biased estimation results (Table 3).

The local spatial autocorrelation can be visually portrayed by plotting the Moran scatter plot. As shown in the Moran scatter plot of PM2.5 in 2014 and 2015 (Figure 2), haze pollution has a significant positive spatial correlation and spatial agglomeration effect. That is, for cities with more serious haze pollution, there are usually one or more cities with more serious haze pollution adjacent to them (i.e., high-high positive correlation). Similarly, cities with lighter haze pollution gather together (i.e., low-low positive correlation). In addition, according to the local Moran’s I, about 110 out of 160 cities have positive local Moran’s I for both PM2.5, further indicating the existence of significant spatial clustering characteristics haze pollution.

Bivariate Moran’s I index is used to explore the spatial correlation between the explained variable and explanatory variables, which is used as a basis for adding spatial lagged terms of explanatory variables in the spatial Durbin model as described earlier (Lee, 2001; Anselin, 2010), the index can be expressed as: I k l = Z k ′ W Z l Z k ′ Z k (12) In Eq. 12, Z _k denotes the normalized Y _k values (explanatory variables) and Z _l denotes the normalized X _l values (explanatory variables), i.e., Z k = [ Y k − μ k ] / σ k , Z l = [ X l − μ l ] / σ l . Y _k and X _l denote the corresponding values, μ _k and μ _l denote the mean values, σ _k and σ _l denote the variance of the corresponding values, and W is the spatial weight matrix normalized by row. Table 4 reports the results.

Before performing the spatial econometric regression, the model is selectively tested by applying the great likelihood estimation and obtaining the likelihood ratio LR. After that, the original hypotheses (H ₀: β = 0 and H ₀: β+αρ = 0) are tested by applying the Wald test and LR test. In addition, the spatial Hausman Test is applied to select fixed effects or random effects. In addition, model selection can also be performed based on the great likelihood value of the model, i.e., Likelihood (Lee and Yu, 2010), and the AIC criterion (Akaike, 1974). Table 5 shows the results, and the Wald test and L.R. test for all three weights significantly reject the original hypotheses of (H ₀: β = 0 and H ₀: β+αρ = 0), indicating that the SDM model is optimal. In addition, the Hausman test for the SDM model with all three weights is significantly positive, rejecting the random effect. Meanwhile, the Likelihood values of the fixed effects are all greater than the random effects, and the AIC values are all smaller than the random effects, so the fixed effects are more suitable. In summary, this paper selected the SDM model with fixed effects for regression analysis.

This paper conducts regression analyses for the POLS model, the spatial random effects SDM model (SDM-RE), and the spatial fixed effects SDM model (SDM-FE) for PM2.5 under three spatial asymmetric weighting schemes, respectively, and the results are shown in Table 6. Among them, according to the identification results based on the Bivariate Moran test, model (1) and model (2) contain the spatial lagged terms W·lnpopdens, W·lnenercons, W·lninvestint, W·lnenerint, W·lnrdint, and W·lntechexp of the explanatory variables. Model (3) and model (4) contain W·lnpgdp, W·lnpopdens, W·lnenercons and W·lnenerint. Model (5) and model (6) contain W·lnenercons, W·lnenerint, and W·lntechexp. Table 6 shows the regression results. Due to the choice of the fixed-effects SDM model in the above model selection, the regression results of model (3), model (5), and model (7) are used in the analysis below.

According to Table 6, the coefficient of W·lnpm passes the significance test at the 10% significance level among the three spatial weighting schemes, which fully indicates that there are significant spatial clustering characteristics and interaction effects of haze pollution in China at present, that is, the increase of PM2.5 concentration in other cities will further aggravate the degree of haze pollution in this city. Under the effect of asymmetric spatial interaction, haze pollution shows a more significant diffusion effect, and the interaction effect of haze pollution among cities becomes increasingly significant. Therefore, combating haze requires regional joint prevention and control by cities, and the effect of haze prevention and control policies of individual cities or a few cities are becoming less and less obvious. Since haze pollution is essentially a socio-economic problem, identifying the socio-economic influences on haze pollution is essential to improving the current haze pollution prevention and control policies.

Based on the regression results below, we further analyze each factor: 1) Economic development. As shown in columns (3), (5), and (7) of Table 6, the primary terms of economic development are all positive, and the secondary terms are all negative. Except for model (5), which is insignificant, model (3) and model (5) pass the 1% significance level, indicating an inverted U-shaped curve between haze pollution and economic development under the asymmetric interaction. It is consistent with the classical EKC hypothesis that the haze pollution level tends to increase and decrease with increasing levels of economic development. On the surface, the haze pollution level decreases after reaching a peak as the economy continues to develop, so there is no need to prevent and control the increasing haze pollution in economic development. However, it is worth noting that we further calculate the inflection point of the inverted U-shaped curve (in the model (7), for example), which is 150.03 million yuan, and only a few of the 160 sample cities exceeded the inflection point. The cities that exceeded the inflection point in 2014 were Daqing, Dongying, Yangquan, Ordos, and Karamay, and in 2015 were Dongying, Daqing, Yantai, Shenzhen, and Ordos, and about 97% of the cities did not cross the inflection point. Therefore, we cannot expect that haze pollution can be improved just by developing the economy, which requires that we should be aware that the management of haze pollution is a long-term and complex process and focus on continuously improving the prevention and control policies of haze pollution while developing the economy. 2) Population density. Observing the regression results of model (3), model (5), and model (7) in Table 6, it can be found that the regression coefficients of population density show significant positive correlations. It indicates that there is a significant positive correlation between population density and haze pollution. In other words, the higher the population density, the more serious the haze pollution. With the accelerated urbanization process, the population is rapidly concentrated in cities, especially in large and medium-sized cities, and the population density is too large to exceed the carrying capacity of the urban environment. In addition, the increase in fossil energy consumption and traffic congestion caused by the over-concentration of the population are also important reasons for the increase in haze pollution. 3) Technological progress. According to Table 6, R&D intensity has a positive impact on haze pollution under all three spatial weighting schemes, but none of them is significant except for the regression coefficient of the model (7). This paper indicates that the current R&D investment in China is not strong enough, and technological progress does not have a significant positive effect on haze pollution. It also indicates that the directional pair of R&D investment is at odds with improving haze pollution. In addition, in terms of science and technology expenditure, the regression coefficients of models (3), (5), and (7) are significantly negative, indicating that the current government financial investment in science and technology has played a beneficial role in improving haze pollution. Therefore, in the future economic development process, the investment in technology dedicated to improving haze pollution should be strengthened, such as increasing the support for new energy and other green technologies and other measures. 4) Energy consumption. The regression coefficients of models (3), (5), and (7) are all positive, and they pass the 1% significance level except for model (7). This result indicates that the current energy consumption structure of China is not reasonable enough, and the coal-based energy consumption structure aggravates the haze pollution. 5) Energy intensity. From the regression results, it can be seen that the regression coefficients are positive and insignificant except for model (5). This indicates that the energy consumption of China’s 10,000 Yuan GDP is gradually decreasing, but it still aggravates haze pollution. Therefore, it is necessary to strengthen energy conservation and emission reduction and enhance the use of clean energy and new energy to further reduce the energy consumption of 10,000 Yuan GDP. 6) Investment intensity. The regression coefficients of investment intensity in models (3), (5), and (7) are significantly negative, indicating that increasing investment in fixed assets is beneficial to improving haze pollution.

Given that the above regression coefficients can only roughly estimate the effect of each factor on haze pollution and cannot accurately explain its direct and indirect effects on haze pollution, this paper, thus, performs effect decomposition on each explanatory variable, as shown in Table 7. Given that the effect decomposition results under the three weight matrices are in the same direction under the same influencing factor, and the spatial weight W₃ can better fit the asymmetric interaction effects of haze pollution, the spatial weight W₃ is used as an example for analysis in the later section.

The effect decomposition results in Table 7 show that: 1) The direct effect of population density (lnpopdens) is 0.425, which is significantly positive, showing that the increase of population density will significantly aggravate the haze pollution in the area. The indirect effect is also positive but insignificant, showing that the positive effect of population concentration in the region on haze pollution in other cities is insignificant, leading to a significant positive total effect. This fully indicates that with the current gradual acceleration of urbanization, the over-saturation of the population in megacities and large cities will exacerbate urban haze pollution while also producing spatial spillover effects on other cities. Although the spatial spillover effect is not too obvious at present, with the further increase of urbanization level, the population of big cities will flow to second-tier and even third-tier cities. Suppose the urban population structure is not reasonably adjusted. In that case, it will also aggravate the haze pollution in other areas. 2) The direct effect of science and technology expenditure (lntechexp) is −0.038, which is significantly negative, while the indirect effect and the total effect are simultaneously negative, indicating that an increase in the proportion of science and technology expenditure in the government’s fiscal expenditure helps to improve the haze pollution in the region and beyond. However, the indirect, and total effects are not significant, which may exist because government fiscal expenditure has certain regional characteristics, and it is difficult to impact other regions. 3) The direct effect, indirect effect, and the total effect of R&D intensity (lnrdint) are significantly positive, but the indirect effect is not significant. 4) The direct and indirect effects of energy consumption (lnenercons) are significantly positive, showing that the unreasonable energy consumption structure will not only aggravate the level of haze pollution in this region but also increase the level of haze pollution in other regions. Hence, the total effect is also significantly positive. 5) The direct effect, indirect effect, and total effect of energy intensity (lnenerint) are all positive but insignificant, indicating that the current energy consumption of 10,000 Yuan GDP in China’s cities is not very reasonable, and further energy saving and emission reduction are needed to reduce the energy consumption of 10,000 Yuan GDP. 6) The direct effect of investment intensity (lninvestint) is significantly negative, which indicates that fixed-asset investment helps to reduce haze pollution in this city. However, its indirect effect is not significant, which shows that the current investment structure is not very reasonable and cannot produce spillover effects.

Examining the above factors together, from the direct effects of each influencing factor, investment intensity and science and technology investment are conducive to reducing haze pollution in the region, while high population density, unreasonable energy consumption structure, and deviation in the direction of R&D intensity exacerbate the level of haze pollution in the region. All factors are insignificant except that energy consumption will significantly aggravate the indirect effects of haze pollution in other cities. Through the above analysis, we can clearly understand that there is still disunity in the prevention and control of haze pollution in China, which is the main reason why the implementation effect of joint prevention and control policies in China is far below the expected level.