Advancing Green TFP Calculation: A Novel Spatiotemporal Econometric Solow Residual Method and Its Application to China’s Urban Industrial Sectors

Abstract

The Solow residual method, traditionally pivotal for calculating total factor productivity (TFP), is typically not applied to green TFP calculations due to its exclusion of undesired outputs. Diverging from traditional approaches and other frontier methodologies such as Data Envelopment Analysis (DEA) and Stochastic Frontier Analysis (SFA), this paper integrates undesired outputs and three types of spatial spillover effects into the conventional Solow framework, thereby creating a new spatiotemporal econometric Solow residual method (STE-SRM). Utilizing this novel method, the study computes the industrial green TFPs for 280 Chinese cities from 2003 to 2019, recalculates these TFPs using DEA-SBM and Bayesian SFA for the same cities and periods, and assesses the accuracy of the STE-SRM-derived TFPs through comparative analysis. Additionally, the paper explores the statistical properties of China’s urban industrial green TFPs as derived from the STE-SRM, employing Dagum’s Gini coefficient and spatial convergence analyses. The findings first indicate that by incorporating undesired outputs and spatial spillover into the Solow residual method, green TFPs are computable in alignment with the traditional Solow logic, although the allocation of per capita inputs and undesired outputs hinges on selecting the optimal empirical production function. Second, China’s urban industrial green TFPs, calculated using the STE-SRM with the spatial Durbin model with mixed effects as the optimal model, show that cities like Huangshan, Fangchenggang, and Sanya have notably higher TFPs, whereas Jincheng, Datong, and Taiyuan display lower TFPs. Third, comparisons of China’s urban industrial green TFP calculations reveal that those derived from the STE-SRM demonstrate broader but more concentrated results, while Bayesian SFA results are narrower and less concentrated, and DEA-SBM findings sit between these extremes. Fourth, the study highlights significant spatial heterogeneity in China’s urban industrial green TFPs across different regions—eastern, central, western, and northeast China—with evident sigma convergence across the urban landscape, though absolute beta convergence is significant only in a limited subset of cities and time periods.

1. Introduction

TFP, or Total Factor Productivity, is a metric that highlights efficiency by subtracting input factors and intermediate inputs from outputs. This calculation offers insights into how effectively resources are being utilized in production processes, especially within specific sectors like industry. By focusing on TFP, analysts can better understand productivity growth, identify areas for improvement, and track changes in efficiency over time. Academic convention typically dictates that TFPs be calculated across two dimensions: the frontier and the non-frontier dimensions [1,2]. The frontier dimension encompasses methods such as Data Envelopment Analysis (DEA) and Stochastic Frontier Analysis (SFA) [3], while the non-frontier dimension includes algebraic or exponential methods [4,5,6] and the Solow residual method [7,8,9,10]. Green Total Factor Productivities (TFPs) are those that incorporate considerations of green inputs or undesirable outputs [11]. The so-called non-desired output refers to the portion of output generated during the production process, such as wastewater, emissions, particulate matter, noise, and other outputs, that do not meet people’s expectations. In essence, the methodologies for calculating green TFPs mirror those of conventional TFPs, with the primary distinction being whether green input factors or undesirable outputs are considered. Currently, green TFPs are predominantly calculated using frontier methods, particularly through DEA, SFA, and their hybrid applications [12,13,14].

The computation of green Total Factor Productivities (TFPs) using Data Envelopment Analysis (DEA) is primarily based on foundational research by Caves et al. (1982) [15] and Färe et al. (1994) [16]. Within the DEA framework, the DEA-Slack Based Measure (DEA-SBM) emerges as a highly favored method for calculating green TFPs. Typically, the green TFPs computed by DEA-SBM are further decomposed into four components: pure efficiency change, pure technology change, scale efficiency change, and technology scale change, using the Malmquist–Luenberger (ML) index. Three significant advancements have been made in the DEA methodology for calculating green TFPs. First, the SBM algorithm has been expanded to include other algorithms such as the BWMRM and BBAM [17]. Second, the ML index has been extended to incorporate additional indices, such as the Malmquist index, the Luenberger index, the Sequential Malmquist–Luenberger index (SML) [18,19,20], and the ML index enhanced with the bootstrap method [21]. Third, the concept of a single frontier has been broadened to encompass multiple frontiers, including general, scale, and regional frontiers [22]. Utilizing these multiple frontiers, DEA has been extended to a three-level meta-frontier DEA framework for the calculation of green TFPs [13,23].

There are numerous shortcomings in calculating green Total Factor Productivities (TFPs) through Data Envelopment Analysis (DEA), including non-proportional changes in desired and undesired outputs, bias from the assumption of constant returns to scale, and the neglect of random disturbances [24]. Stochastic Frontier Analysis (SFA) can mitigate these issues to a certain extent. The methods within the SFA framework, largely developed and advanced based on the work of Battese et al. (1995) [25], typically break down green TFPs into several components, including the rate of frontier technological progress, the rate of change in technology inefficiency, the rate of change in factor returns to scale, and the rate of change in factor allocation efficiency [26]. Recent advancements in SFA methodology primarily focus on the detailed breakdown of core explanatory variables within the frontier production function, particularly the disaggregation of inputs such as capital and labor. For instance, Zhu et al. (2020) categorized capital into three types: construction and installation engineering, equipment and tools, and other investment types [12]. They also classified labor into skilled and unskilled categories, enabling a more nuanced analysis of input contributions to productivity.

Although Stochastic Frontier Analysis (SFA) addresses some of the limitations of Data Envelopment Analysis (DEA), it has its own shortcomings. These include a relatively limited variety of stochastic frontier function models and the potential for errors in model specification. Recognizing the limitations inherent in both DEA and SFA methodologies, researchers have begun to develop hybrid approaches that integrate both techniques, such as the three-stage DEA analysis [14,27]. This three-stage DEA analysis differs from the previously mentioned three-levels meta-frontier DEA. While the latter expands from a single frontier to multiple frontiers, the three-stage DEA focuses on extending the analysis from a single step to multiple steps, involving an initial DEA, followed by SFA, and then DEA again. In this three-stage approach, the first step calculates the green TFPs and their input slackness using DEA. The second step removes the influences of management inefficiency factors, environmental factors, and random errors from the input slackness, and reconstructs the adjusted inputs based on SFA. The third step then recalculates the green TFPs using these adjusted inputs and another round of DEA analysis. This sequential method aims to enhance the precision and reliability of green TFP measurements by integrating the strengths of both DEA and SFA.

There are limited methods for calculating green Total Factor Productivities (TFPs) from the non-frontier perspective, particularly through the Solow residual method. Recent studies, such as those by Cheng and Jin (2022) and Song et al. (2022) [28,29], have extensively examined the correlations between green TFPs and factors like foreign direct investment, economic agglomeration, and climate change, often employing the Cobb–Douglas production function for analysis. In these studies, green TFPs are typically defined as the output surplus remaining after accounting for inputs and their shares, and are influenced by factors such as foreign direct investment and climate change [30]. However, a critical discrepancy arises in these analyses: the theoretical definition of green TFPs as derived from the conceptual framework of the Solow residual method is often disconnected from their actual empirical measurement. In practice, green TFPs are calculated using methods like the directional distance function (DDF), DEA, or SFA, leading to a significant “two skins” problem where theoretical and practical approaches do not align [29]. This disconnection underscores an urgent need for innovation in the methodology for calculating green TFPs from the Solow residual perspective. Furthermore, the prevailing misconceptions in current research regarding the relationship between green TFPs and other economic or social development variables also highlight the critical necessity for innovative and advanced research that adheres to the core principles of the Solow residual method.

This paper seeks to enhance the methodology for calculating green Total Factor Productivities (TFPs) by integrating undesired outputs and spatiotemporal econometric models into the traditional Solow residual method, thereby creating a novel spatiotemporal econometric Solow residual method (STE-SRM). This new approach is applied to compute the industrial green TFPs of 280 major Chinese cities from 2003 to 2019. Additionally, the accuracy of this novel method will be evaluated through comparative analysis with results obtained from the DEA-SBM and Bayesian SFA. The statistical characteristics of the industrial green TFPs calculated by the new method will also be explored, utilizing Dagum’s Gini coefficient and convergence analyses, including sigma and beta convergence. This innovation is expected to significantly enhance the application of the Solow residual method in calculating industrial green TFPs and contribute to advancing the methodology for green TFP calculation from the non-frontier perspective.

2. Redefining Green TFP Calculation: Innovations from the Solow Residual Method to the Spatiotemporal Econometric Solow Residual Method (STE-SRM)

2.1. Preliminary Extension of the Solow Residual Method: Incorporating Undesired Outputs

The Solow residual method is a widely utilized technique for Total Factor Productivity (TFP) calculation. Its fundamental principle involves determining the ratio of per capita outputs to per capita inputs and their respective shares. Operating under assumptions of neutral technological progress and constant returns to scale, this method typically employs the Cobb–Douglas production function, represented as $Y=A{K}^{\alpha }{L}^{\beta }$, and the Solow residual method calculates Total Factor Productivity (TFP) using the formula $A=y/k^{\alpha }$. In the production function, Y, A, K, and L represent total output, technology, capital, and labor, respectively. Conversely, y and k are the per capita measures of output and capital. Additionally, $\alpha$ signifies the input share of capital.

Green Total Factor Productivities (TFPs) account for undesired outputs. By incorporating these outputs as environmental inputs within the Cobb–Douglas production function, we can calculate green TFPs in a manner akin to previous analyses. This adaptation involves redefining the Cobb–Douglas production function to include undesired outputs, as shown in Equation (1).

$$\left({\displaystyle \prod _{j=1}^J{Y}_j^{{\beta }_j}}\right)Y^{{\beta }_0}=A{L}^{{\alpha }_0}{\displaystyle \prod _{q=1}^Q{X}_q^{{\alpha }_q}}$$

(1)

In Equation (1), A represents neutral technological progress; L denotes labor, and ${\alpha }_0$ represents the input share of labor; $X_q$ represents the $q^{th}$ inputs, $q=1,2,\cdots ,Q$, ${\alpha }_q$ denote the share of $X_q$, and ${\displaystyle \prod \left(·\right)}$ is a symbol representing a continuous multiplicative factor. Under the assumption of constant returns to scale, ${\alpha }_0+{\displaystyle \sum _{q=1}^Q{\alpha }_q}=1$. In Equation (1), Y indicates desired outputs, while $Y_j$ indicates the $j^{th}$ undesired outputs, $j=1,2,\cdots ,J$; ${\beta }_0$ and ${\beta }_j$ indicate the share of desired outputs and the $j^{th}$ undesired outputs, respectively. We define ${\beta }_0+{\displaystyle \sum _{j=1}^J{\beta }_j}=1$, then ${\beta }_0$ presents the share of the desired outputs in the constitution of the whole outputs. A higher value of ${\beta }_0$ means better efficiency of outputs. Equation (1) can be transformed into Equation (2) if the absolute inputs and outputs are replaced with per capita values.

$$\left[{\displaystyle \prod _{j=1}^J{\left(\frac{Y_j}{Y}\right)}^{{\beta }_j}}\right]\left(\frac{Y}{L}\right)=A{\displaystyle \prod _{q=1}^Q{\left(\frac{X_q}{L}\right)}^{{\alpha }_q}}$$

(2)

We define $yy=Y/L$, $x{x}_q=X_q/L$, $y{y}_j=Y_j/Y$, the production function with Log form can be shown as in Equation (3), where ${\epsilon }_1$ is the random disturbance term,${\epsilon }_1\sim I.I.D\left(0,{\sigma }_1^2\right)$.

$$Ln\left(yy\right)=Ln\left(A\right)+{\displaystyle \sum _{q=1}^Q{\alpha }_{q}Ln\left(x{x}_q\right)}-{\displaystyle \sum _{j=1}^J{\beta }_{j}Ln\left(y{y}_j\right)}+{\epsilon }_1$$

(3)

From the parameters estimated in Equation (3), the green TFP can be determined as outlined in Equation (4), where $GTF{P}_0$ represents the green TFP that incorporates undesired outputs. In Equation (4), the symbol of ^ indicates the estimators of the corresponding parameters, while other symbols are defined the same as in Equation (3).

$$GTF{P}_0=\frac{yy}{{\displaystyle \prod _{q=1}^Q{\left(x{x}_q\right)}^{{\widehat{\alpha }}_q}}{\displaystyle \prod _{j=1}^J{\left(y{y}_j\right)}^{{\widehat{\beta }}_j}}}$$

(4)

2.2. Further Extension: Defining the STE-SRM with Consideration of Three Types of Spillover Effects

Further discussion is warranted if the spillover effects from spatial neighbors are integrated into the production function. The general nesting spatial model, a comprehensive formulation in spatial econometrics, encapsulates various types of spillover effects emanating from the outputs, inputs, or disturbance factors of neighboring regions [31,32]. By embedding this model into the logarithmic form of the production function in Equation (3), it can be transformed into a revised version, as detailed in Equations (5) and (6).

$$y=\rho STWy+X\beta +STWX\theta +u+v+{\mu }_1$$

(5)

$${\mu }_1=\lambda \left(STW\times {\mu }_1\right)+{\epsilon }_2$$

(6)

In Equations (5) and (6), the parameters are set such that $y=Ln\left(yy\right)$ and $X=\left[\begin{array}{cccccc}Ln\left(x{x}_1\right)& \cdots & Ln\left(x{x}_Q\right)& Ln\left(y{y}_1\right)& \cdots & Ln\left(y{y}_J\right)\end{array}\right]$; the spatial-temporal weight matrix (STW) is typically constructed as the Kronecker product of the temporal weight matrix and the spatial weight matrix, effectively integrating both spatial and temporal dimensions; in the model, $\rho$ and $\lambda$ represent spatial correlation coefficients; $\beta$ and $\theta$ are exogenous parameters of the independent variables, $\beta ={\left[\begin{array}{cccccc}{\alpha }_1& \cdots & {\alpha }_Q& -{\beta }_1& \cdots & -{\beta }_J\end{array}\right]}^{\prime }$, $\theta ={\left[\begin{array}{cccccc}{\theta }_1& \cdots & {\theta }_Q& -{\gamma }_1& \cdots & -{\gamma }_J\end{array}\right]}^{\prime }$; $u$ and $v$ represent the period effect and individual effect, respectively. These effects can either be fixed or random, depending on the specific requirements and assumptions of the analysis; ${\mu }_1$ and ${\epsilon }_2$ are disturbance terms, ${\epsilon }_2\sim I.I.D\left(0,{\sigma }_2^2\right)$, and the distribution of ${\mu }_1$ is decided by Equation (6).

Based on the framework outlined in Equations (5) and (6), and considering various hypotheses such as $\rho =0$, or $\lambda =0$, or $\theta =0$, a total of seven reduced models can be derived from the original specifications. These seven models are typically referred to as the non-spatial model (NSM), the spatial X-lag model (SXL), the spatial autoregressive model (SAR), the spatial error model (SEM), the spatial Durbin model (SDM), the spatial autocorrelation model (SAC), and the spatial Durbin error model (SDEM), as described by Zhao et al. (2022) [33]. In the aforementioned eight models, the parameters, the variance of the disturbance terms, and the variance-covariance matrices can be estimated using the robust analytical framework of spatial econometrics pioneered by LeSage (2009) and Elhorst (2014) [31,32]. This estimation methodology integrates maximum likelihood estimation, the Generalized Method of Moments, Bayesian Markov Chain Monte Carlo simulation techniques, and other approaches.

The green Total Factor Productivity (TFP), which accounts for both undesired outputs and spillover effects, can be calculated using Equation (7). This equation utilizes the data-generating process from the optimal model selected among the estimated models from Equations (5) and (6) and their corresponding reduced models.

$$GTFP=\frac{yy}{{\displaystyle \prod _{q=1}^Q{\left(x{x}_q\right)}^{{\widehat{\eta }}_q}}{\displaystyle \prod _{j=1}^J{\left(y{y}_j\right)}^{{\widehat{\tau }}_j}}}$$

(7)

In Equation (7), $GTFP$ is the green TFP considering both the undesired outputs and spillover effects; the definition of $yy$, $x{x}_q$, and $y{y}_j$ is the same as in Equation (4); ${\widehat{\eta }}_q$ is the estimated share of the $q^{th}$ inputs per capita, ${\widehat{\eta }}_q=\frac{1}{NT}trace\left[\widehat{\Theta }\left({\widehat{\alpha }}_q{I}_{NT}+{\widehat{\theta }}_{q}STW\right)\right]$, where $\widehat{\Theta }={\left(I_{NT}-\widehat{\rho }STW\right)}^{-1}$, N and T respectively indicate the total number of the spatial units and periods, $I_{NT}$ is the unit matrix with both N*T rows and columns, and $trace\left(·\right)$ is the trace statistic; ${\widehat{\tau }}_j$ is the estimated share of the $j^{th}$ undesired outputs in the constitution of the whole outputs, ${\widehat{\tau }}_j=\frac{1}{NT}trace\left[-\widehat{\Theta }\left({\widehat{\beta }}_j{I}_{NT}+{\widehat{\gamma }}_{j}STW\right)\right]$; the symbol of ^ still indicates the estimators.

The green Total Factor Productivity (TFP) calculated using Equation (7) resembles that derived from Equation (4). However, the TFP from Equation (7) is deemed more significant and valuable due to its consideration of the spatiotemporal impacts from spatial neighbors. If the true values of $\rho$, $\lambda$, and $\theta$ are zero, the green TFPs calculated from both equations would be equivalent. Given this context, the method of calculating green TFP using Equations (5) and (6) can be termed the Spatiotemporal Econometric Solow Residual Method (STE-SRM). This designation is appropriate because it adheres to the core analytical logic of the Solow residual method while incorporating global spatiotemporal econometric models into the analysis.

3. Applying the New Method to Calculate China’s Urban Industrial Green TFPs

Industry is a crucial component of the modern industrial system and plays a pivotal role as both a carrier and cornerstone in advancing the construction of a manufacturing powerhouse. In recent years, the share of industrial-added value in China’s GDP has been on a decline. At the beginning of 2024, this proportion stood at 33.19%, which is 11.67% lower than in 2003 when China launched its new industrialization strategy. Amidst the ongoing reduction in the relative scale of industrial development, fully enhancing industrial efficiency and fostering the green transformation of the sector are perpetual themes in the high-quality development of China’s economy. Chinese cities serve as vital spatial units for economic and particularly industrial development. The industrial green Total Factor Productivities (TFPs) in these cities are integral to the broader national economic system and the pursuit of sustainable industrial growth. This section will detail the calculation of China’s urban industrial green TFPs using the innovative Spatiotemporal Econometric Solow Residual Method (STE-SRM).

3.1. Configuration of the Empirical Production Function Model

The STE-SRM offers enhanced accuracy because it incorporates the spatiotemporal spillover effects from various spatial units into the calculation of green TFPs. In this section, the STE-SRM will be utilized to calculate the green TFPs of China’s urban industrial sectors. In the production model for China’s urban industries, three types of inputs and three types of undesired outputs are considered. The input factors include capital ($K$), labor ($L$), and energy ($E$), while the undesired outputs comprise industrial wastewater emissions ($Y_{WW}$), industrial sulfur dioxide emissions ($Y_{S{O}_2}$), and industrial dust emissions ($Y_{S\&D}$). We define $kk=K/L$, $ee=E/L$, $y_{ww}=Y_{WW}/Y$, $y_{s{o}_2}=Y_{S{O}_2}/Y$, $y_{s\&d}=Y_{S\&D}/Y$, where the empirical production function model for China’s urban industrial sectors is specified in Equations (8) and (9).

$$\begin{array}{c}Ln\left(yy\right)=Ln\left(A\right)+\rho \left[STW\times Ln\left(yy\right)\right]+{\alpha }_{1}Ln\left(kk\right)+{\alpha }_{2}Ln\left(ee\right)-{\beta }_{1}Ln\left(y_{ww}\right)-{\beta }_{2}Ln\left(y_{s{o}_2}\right)\\ -{\beta }_{3}Ln\left(y_{s\&d}\right)+{\theta }_1\left[STW\times Ln\left(kk\right)\right]+{\theta }_2\left[STW\times Ln\left(ee\right)\right]-{\gamma }_1\left[STW\times Ln\left(y_{ww}\right)\right]\\ -{\gamma }_2\left[STW\times Ln\left(y_{s{o}_2}\right)\right]-{\gamma }_3\left[STW\times Ln\left(y_{s\&d}\right)\right]+u+v+{\mu }_2\end{array}$$

(8)

$${\mu }_2=\lambda STW{\mu }_2+{\epsilon }_3$$

(9)

In Equations (8) and (9), $\rho$ and $\lambda$ are spatial correlation coefficients; $\alpha$, $\beta$, $\theta$ and $\gamma$ are exogenous parameters similar to those in Equations (5) and (6), and Q and J will be 2 and 3, respectively; $u$ and $v$ are the period effect and individual effect, respectively; ${\mu }_2$ and ${\epsilon }_3$ are disturbance terms, ${\epsilon }_3~I.I.D\left(0,{\sigma }_3^2\right)$, and the distribution of ${\mu }_2$ is decided by Equation (9); $STW$ is the spatiotemporal weight matrix that represents the spatial and temporal spillover effects of the spatial units during the research periods. The construction of this matrix will be elaborated on in a subsequent section.

3.2. Variables and Data Overview

This paper focuses on 280 cities in mainland China as the spatial units and covers the period from 2003 to 2019. The spatial units include 15 sub-provincial cities: Guangzhou, Wuhan, Harbin, Shenyang, Chengdu, Nanjing, Xi’an, Changchun, Jinan, Hangzhou, Dalian, Qingdao, Shenzhen, Xiamen, and Ningbo. Additionally, the study encompasses all the prefectural cities in mainland China, with the exception of Danzhou and Sansha in Hainan, Bijie and Tongren in Guizhou, Haidong in Qinghai, Turpan and Hami in Xinjiang, and Xigaze, Changdu, Linzhi, Shannan, and Naqu in Tibet. These 12 cities are excluded due to adjustments in their administrative levels between 2003 and 2019. The selection of the research period is based on two factors: 2003 marks a significant milestone in the initiation of China’s new industrialization strategy, and data post-2020 has not yet been updated.

This paper utilizes various significant statistical indices to illustrate the inputs and outputs of China’s urban industrial sectors. Specifically, the gross value of industrial outputs (in CNY 100 million) is used to represent industrial outputs. Inputs are quantified using three measures: industrial capital stock (in CNY 100 million), calculated via the perpetual inventory method; total industrial employment (in ten thousand people); and industrial electricity consumption (in 100 million kWh). Additionally, three types of undesired outputs are tracked: industrial wastewater emissions (in ten thousand tons), industrial sulfur dioxide emissions (in ten thousand tons), and smoke (dust) emissions (in ten thousand tons). It is important to note that all data concerning inputs and outputs are collected at the city level for the cities mentioned previously. Most of the data is accessible from several sources, including the EPS statistical database, China Economic and Social Development Statistical Database, DRCnet Statistical Database, and published annual reports such as the China Urban Statistical Yearbook, China Regional Statistical Yearbook, and the China Statistical Yearbook.

In this paper, two important techniques are employed for data processing. First, the gross value of industrial outputs and industrial capital stock are converted to real values (constant 1990 prices) using the GDP deflator and the fixed assets investment price index, respectively. Second, missing data are interpolated using the method of a five-year moving average. Additionally, the interpolated data are adjusted based on the development trends of neighboring cities within the same provinces. Further details on the data processing methods are thoroughly documented in previously published papers by the author. The descriptive statistical characteristics of the data are presented in

Table 1. Descriptive statistics of related variables (logarithmic form) for calculating China’s urban industrial green TFPs.

	$\mathit{L}\mathit{n}\left(\mathit{y}\mathit{y}\right)$	$\mathit{L}\mathit{n}\left({\mathit{y}}_{\mathit{w}\mathit{w}}\right)$	$\mathit{L}\mathit{n}\left({\mathit{y}}_{\mathit{s}{\mathit{o}}_2}\right)$	$\mathit{L}\mathit{n}\left({\mathit{y}}_{\mathit{s}\&\mathit{d}}\right)$	$\mathit{L}\mathit{n}\left(\mathit{k}\mathit{k}\right)$	$\mathit{L}\mathit{n}\left(\mathit{e}\mathit{e}\right)$
Mean	3.5893	2.4134	-4.7971	−5.1456	2.6460	1.8661
Median	3.6953	2.3427	−4.7338	−5.1417	2.6632	1.8044
Maximum	6.2935	6.9681	−0.4706	3.2454	5.9188	5.9422
Minimum	0.1778	−0.8902	−11.6014	−12.3200	0.2191	−3.5941
Std. Dev.	0.7914	1.1041	1.4821	1.6388	0.8355	0.9109
Skewness	−0.4191	0.3424	−0.3999	0.0775	0.1042	0.4612
Kurtosis	3.2116	3.0585	3.4107	3.8777	2.9826	4.8804
Jarque-Bera	148.2233	93.6619	160.3404	157.5696	8.6778	870.0061
Probability	0.0000	0.0000	0.0000	0.0000	0.0131	0.0000
Sum	17,084.96	11,487.66	−22,834.20	−24,492.83	12,595.08	8882.55
Sum Sq. Dev.	2980.62	5801.75	10,453.91	12,781.42	3322.28	3948.82
Observations	4760	4760	4760	4760	4760	4760

Note: the results were obtained based on MATLAB R2023a and EViews 11.0.

.

3.3. Development of the Spatiotemporal Weight Matrix

The spatiotemporal weight matrix ($STW$) is crucial for estimating the empirical production function as outlined in Equations (8) and (9). This matrix is typically constructed through the Kronecker product of the temporal weight matrix ($TW$) and the spatial weight matrix ($W$), denoted as $STW=Kron\left(TW,W\right)$, where $Kron\left(·\right)$ represents the Kronecker product operation. The spatial weight matrix reflects the spillover effects among spatial units in a special year, and it is also the foundation for the construction of the temporal weight matrix. Thus, construction of the spatiotemporal weight matrix will begin with the designs of the elements of $W$.

Since the 280 cities included in this study do not form a complete set, the spatial weight matrix $W$ is not designed based on spatial proximity alone. Instead, it is determined by the reciprocal of the square of the distance in longitude and latitude between pairs of cities. Subsequently, the elements of $W$ are adjusted according to the boundaries of cities that have effective impacts. According to research by Guo and Fan (2022) [34], the adaptive bandwidth for China’s urban industrial production is 46. This means that elements in $W$ are set to zero if they are ranked lower than 46 in each row, and are retained if their ranks are within the top 46 in each row. Following these adjustments, $W$ is row-stochastic standardized to ensure that the sum of each row equals one, as depicted in Figure 1a. It is important to note that the order of the 280 cities significantly influences the structure of the spatial weight matrix. For this study, cities are ranked according to their positions within each province as listed in the EPS statistical database. The order of the provinces is as follows: Hebei, Shanxi, Inner Mongolia, Liaoning, Jilin, Heilongjiang, Jiangsu, Zhejiang, Anhui, Fujian, Jiangxi, Shandong, Henan, Hubei, Hunan, Guangdong, Guangxi, Hainan, Sichuan, Guizhou, Yunnan, Shaanxi, Gansu, Qinghai, Ningxia, and Xinjiang.

The temporal weight matrix captures the variations in spillover effects among spatial units across different years, and this study employs the ratios of global Moran’s I to illustrate these changes. Typically, the global Moran’s I index is derived from the standardized spatial weight matrix and the data for the dependent variable or one of the independent variables. However, this method of index calculation can introduce endogeneity, which may affect the validity of the results. To avoid possible endogeneity in the designs of the elements of $TW$, this paper calculates the global Moran’s I using the estimated residuals (${\widehat{\epsilon }}_4$) in the following Equation (10), as shown in Figure 1b. The definitions of the variables and parameters in Equation (10) are the same as in Equations (8) and (9), and ${\epsilon }_4$ is also the random disturbance term, ${\epsilon }_4~IID\left(0,{\sigma }_4^2\right)$.

$$Ln\left(yy\right)=Ln\left(A\right)+{\alpha }_{1}Ln\left(kk\right)+{\alpha }_{2}Ln\left(ee\right)-{\beta }_{1}Ln\left(y_{ww}\right)-{\beta }_{2}Ln\left(y_{s{o}_2}\right)-{\beta }_{3}Ln\left(y_{s\&d}\right)+{\epsilon }_4$$

(10)

Based on the global Moran’s I indices depicted in Figure 1b, the temporal weight matrix is determined, as illustrated in Figure 1c. From this figure, it is evident that the temporal weight matrix is a lower triangular matrix with the number of rows and columns both equal to T, where T represents the total number of periods analyzed. The characteristics of this lower triangular structure are due to the fact that temporal spillover effects originating from a specific year typically exhibit a time lag. Utilizing the Kronecker product of the temporal weight matrix $TW$ from Figure 1c and $W$ from Figure 1a, the spatiotemporal weight matrix is constructed, as shown in Figure 1d. This comprehensive matrix integrates both spatial and temporal dynamics, enabling a nuanced analysis of spillover effects across different times and locations.

3.4. Estimation and Model Selection for the Empirical Production Function

Using the spatiotemporal weight matrix shown in Figure 1d and integrating the variables and data described in Section 3.2, the models from Equations (8) and (9), along with their degraded versions, can be estimated. The results, as presented in

Table 2. Estimated results of eight spatial models without fixed or random effects.

	NSM	SXL	SAR	SEM	SDM	SDEM	SAC	GNSM
Const.	1.4267 (27.1) ***	0.8689 (7.91) ***	0.7545 (12.71) ***	2.3353 (20.04) ***	0.3893 (3.99) ***	2.9512 (11.30) ***	4.5441 (80.90) ***	3.3353 (12.12) ***
$Ln\left(kk\right)$	0.3254 (31.53) ***	0.3864 (35.13) ***	0.2919 (29.27) ***	0.3612 (36.96) ***	0.3830 (39.49) ***	0.3699 (37.43) ***	0.3812 (40.04) ***	0.3691 (37.29) ***
$Ln\left(ee\right)$	0.2232 (26.02) ***	0.2200 (24.87) ***	0.1887 (22.61) ***	0.1920 (24.16) ***	0.1991 (25.46) ***	0.2090 (25.43) ***	0.1928 (24.73) ***	0.2107 (25.66) ***
$Ln\left(y_{ww}\right)$	−0.1042 (−12.64) ***	−0.1580 (−17.40) ***	−0.0925 (−11.73) ***	−0.1402 (−16.79) ***	−0.1456 (−18.17) ***	−0.1474 (−17.49) ***	−0.1516 (−18.78) ***	−0.1486 (−18.26) ***
$Ln\left(y_{s{o}_2}\right)$	−0.0701 (−8.87) ***	−0.0519 (−6.29) ***	−0.0340 (−4.40) ***	−0.0260 (−3.59) ***	−0.0413 (−5.67) ***	−0.0296 (−3.99) ***	−0.0341 (−4.89) ***	−0.0297 (−4.05) ***
$Ln\left(y_{s\&d}\right)$	−0.1555 (−23.73) ***	−0.1306 (−19.27) ***	−0.1138 (−17.37) ***	−0.1008 (−16.52) ***	−0.1130 (−18.85) ***	−0.1211 (−19.53) ***	−0.1057 (−17.80) ***	−0.1242 (−20.10) ***
$STW\times Ln\left(kk\right)$		−0.2647 (−10.16) ***			−0.5436 (−22.56) ***	−0.2732 (−8.01) ***		−0.2158 (−4.98) ***
$STW\times Ln\left(ee\right)$		0.1106 (5.36) ***			−0.1152 (−6.02) ***	0.1230 (4.93) ***		0.1544 (5.65) ***
$STW\times Ln\left(y_{ww}\right)$		0.1696 (9.09) ***			0.1881 (11.43) ***	0.0784 (2.51) **		0.0398 (1.38)
$STW\times Ln\left(y_{s{o}_2}\right)$		0.0369 (1.54)			0.1702 (7.96) ***	0.2359 (9.77) ***		0.2204 (9.09) ***
$STW\times Ln\left(y_{s\&d}\right)$		−0.1773 (−11.16) ***			−0.0097 (−0.66)	−0.1866 (−8.36) ***		−0.2097 (−9.50) ***
$\rho$			0.3780 (21.27) ***		0.9010 (38.57) ***		−0.3050 (−89.77) ***	−0.2000 (−2.80) ***
$\lambda$				0.8860 (49.16) ***		0.9120 (85.53) ***	0.9400 (184.88) ***	0.9000 (74.52) ***
${\widehat{R}}^2$	0.7101	0.7308	0.6928	0.7785	−0.4350	0.7909	0.7849	0.7918
${\sigma }^2$	0.1816	0.1686	0.1653	0.1386	0.1311	0.1307	0.1346	0.1301
$Log\left(L\right)$	−2690.40	−2511.90	−830.05	−430.42	−328.35	−293.62	−403.77	−1932.70

Note: the above outputs were collected based on MATLAB R2023a. ( ) denotes the T-statistic, *** and ** indicate that the hypothesis tests were passed at significance levels of 1% and 5%, respectively.

, indicate the absence of individual or period effects, suggesting that $u=0$ and $v=0$. From the analysis presented in Table 2, it is evident that the Spatial X-lag Model (SXL) and the General Nested Spatial Model (GNSM) are not optimal choices. This conclusion is based on the inability of some estimated parameters in these models to pass hypothesis tests and achieve statistical significance. Similarly, the Spatial Durbin Model (SDM) is not deemed optimal due to its negative goodness of fit and the lack of significance in some of its parameters. Furthermore, although all parameters of the Spatial Autocorrelation Model (SAC) are statistically significant, this model is also not suitable because its spatial correlation coefficient $\rho$ is negative, which contradicts the economic implications suggested by the positive global Moran’s I indices shown in Figure 1b. Likewise, the Non-Spatial Model (NSM) is not considered optimal as it lacks spatial spillover terms, despite the significant positive values of the global Moran’s I indices, indicating the importance of spatial effects in the data.

Three models remain under consideration: the Spatial Autoregressive Model (SAR), the Spatial Error Model (SEM), and the Spatial Durbin Error Model (SDEM). This paper employs Lagrange Multiplier (LM) tests and Likelihood Ratio (LR) tests to evaluate these models, assessing their statistical robustness and suitability for the analysis presented. For the model selection between SAR and SEM, LM tests play an important role, wherein $LM\_Lag$ and $LM\_Error$ are 840.42 and 2576.83, respectively, bigger than the threshold of the Chi-squared test at the level of 1% (6.64), while $Robust\_LM\_Lag$ and $Robust\_LM\_Error$ are 1.13 and 1737.53, respectively. The LM tests show that the SEM will be better than the SAR. Further selection can be carried out between SEM and SDEM by LR tests. Under the null hypothesis of $H_0:{\theta }_1={\theta }_2={\gamma }_1={\gamma }_2={\gamma }_3$, the LR statistic is 273.6, bigger than the threshold of the Chi-squared test with a freedom of 5 and a significance level at 1%. This suggests that the Spatial Durbin Error Model (SDEM) may perform better than the Spatial Error Model (SEM) and could potentially be the optimal model for this analysis.

Further analyses to select among the Spatial Autoregressive Model (SAR), Spatial Error Model (SEM), and Spatial Durbin Error Model (SDEM) include consideration of these models with either individual or period effects. Given that the majority of cities are included in the analysis of China’s urban industrial sectors, it is deemed unnecessary to focus solely on models with random effects [35]. Therefore, versions of the SAR, SEM, and SDEM that incorporate fixed effects—whether individual, period, or both—are estimated, as detailed in

Table 3. Estimated results of SAR, SEM, and SDEM Models with fixed effects.

	SAR			SEM			SDEM
	Individual Fixed	Period Fixed	Both Fixed	Individual Fixed	Period Fixed	Both Fixed	Individual Fixed	Period Fixed	Both Fixed
Const.	0.0775 (21.17) ***	−0.0117 (−2.19) **	−1.3077 (−24.98) ***	0.2537 (7.90) ***	−0.2291 (−4.63) ***	−1.6788 (−100.74) ***	0.0643 (1.80) *	−0.2107 (−4.18) ***	−1.6166 (−74.71) ***
$Ln\left(kk\right)$	0.4159 (39.18) ***	0.3643 (37.74) ***	0.5454 (53.53) ***	0.4329 (41.42) ***	0.4108 (43.12) ***	0.5508 (77.17) ***	0.4841 (44.26) ***	0.4057 (41.74) ***	0.5468 (54.21) ***
$Ln\left(ee\right)$	0.1918 (20.29) ***	0.1924 (24.98) ***	0.1345 (16.99) ***	0.1643 (18.14) ***	0.1943 (26.01) ***	0.1313 (16.65) ***	0.1712 (19.98) ***	0.2054 (25.89) ***	0.1326 (16.80) ***
$Ln\left(y_{ww}\right)$	−0.1986 (−23.47) ***	−0.1036 (−14.06) ***	−0.1725 (−24.17) ***	−0.1788 (−21.90) ***	−0.1437 (−18.67) ***	−0.1722 (−30.79) ***	−0.1907 (−24.29) ***	−0.1390 (−17.86) ***	−0.1720 (−24.93) ***
$Ln\left(y_{s{o}_2}\right)$	0.0633 (9.08) ***	−0.1139 (−15.28) ***	−0.0670 (−10.21) ***	0.0432 (6.42) ***	−0.0886 (−12.50) ***	−0.0665 (−11.01) ***	−0.0100 (−1.45)	−0.0844 (−11.33) ***	−0.0667 (−10.28) ***
$Ln\left(y_{s\&d}\right)$	−0.0932 (−14.43) ***	−0.0869 (−14.49) ***	−0.0581 (−10.71) ***	−0.0816 (−13.18) ***	−0.0902 (−15.87) ***	−0.0596 (−11.44) ***	−0.0851 (−14.43) ***	−0.0993 (−16.57) ***	−0.0589 (−10.91) ***
$STW\times Ln\left(kk\right)$							−0.4753 (−12.45) ***	−0.0360 (−1.04)	−0.0318 (−1.00)
$STW\times Ln\left(ee\right)$							0.3861 (10.46) ***	0.1056 (4.39) ***	0.1172 (3.54) ***
$STW\times Ln\left(y_{ww}\right)$							0.0074 (0.23)	0.0020 (0.07)	−0.0367 (−2.16) **
$STW\times Ln\left(y_{s{o}_2}\right)$							0.3875 (13.05) ***	0.0404 (1.67) *	0.0108 (0.47)
$STW\times Ln\left(y_{s\&d}\right)$							−0.2241 (−7.90) ***	−0.1131 (−5.35) ***	−0.0152 (−0.61)
$\rho$	0.1620 (66.49) ***	0.3850 (22.94) ***	0.1070 (53.95) ***
$\lambda$				0.7230 (35.54) ***	0.8730 (50.53) ***	0.3310 (87.89) ***	0.7520 (42.41) ***	0.8740 (43.91) ***	0.3100 (45.23)
${\widehat{R}}^2$	0.8032	0.5725	0.6043	0.8277	0.6805	0.6096	0.8464	0.6836	0.6105
${\sigma }^2$	0.0702	0.1339	0.0449	0.0639	0.1180	0.0446	0.0569	0.1167	0.0444
$Log\left(L\right)$	1217.00	−328.64	2281.40	1426.10	−45.60	2295.00	1700.40	−19.45	2303.40

Note: the above outputs were collected based on MATLAB R2023a. ( ) denotes the T-statistic, ***, **, and * mean having passed the hypothesis test with a significance level of 1%, 5%, and 10%, respectively.

. Models incorporating individual fixed effects and those with both fixed effects show positive $Log\left(L\right)$ values in Table 3, indicating that these configurations may not be optimal. For the models with period fixed effects, the goodness of fit for SAR and SEM is lower than their counterparts in Table 2, and some parameters in the SDEM with period fixed effects lack significance. Thus, compared to their versions without fixed effects, the models with fixed effects do not exhibit improved estimation results. In conclusion, the SDEM without fixed or random effects emerges as the optimal model for analyzing China’s urban industrial production. This model will also be used later to calculate China’s urban industrial green TFPs.

3.5. Calculation of China’s Urban Industrial Green TFPs Using the STE-SRM

Based on Equation (7) and the estimates derived from the Spatial Durbin Error Model (SDEM) presented in Table 2, it is possible to compute the green total factor productivity (TFP) of China’s urban industries, as depicted in Figure 2. It is important to highlight that the SDEM, functioning as a degradation model in accordance with Equations (8) and (9), operates under the condition where $\widehat{\rho }=0$, and then $\widehat{\Theta }=I_{NT}$; thus, the estimated share of inputs per capita in Equation (7) will be changed to ${\widehat{\eta }}_q=\frac{1}{NT}trace\left({\widehat{\alpha }}_q{I}_{NT}+{\widehat{\theta }}_{q}STW\right)$, while the estimated share of the $j^{th}$ undesired outputs in the constitution of the whole outputs will be also changed to ${\widehat{\tau }}_j=\frac{1}{NT}trace\left[-\left({\widehat{\beta }}_j{I}_{NT}+{\widehat{\gamma }}_{j}STW\right)\right]$ correspondingly. The input shares of per capita industrial capital and per capita industrial energy in China are 0.3699 and 0.2090, respectively. This distribution suggests that China’s industrial growth remains significantly dependent on labor inputs. Additionally, among the three undesired outputs—industrial wastewater emissions, industrial sulfur dioxide emissions, and industrial smoke (dust) emissions—their respective shares are 0.1474, 0.0296, and 0.1211. Consequently, these figures imply that the proportion of expected output in China’s urban industrial sectors is approximately 0.7019.

Furthermore, this study identifies the 10th and 90th percentiles of all industrial Green TFPs from 2003 to 2019 as critical benchmarks and details the cities whose industrial green TFPs are in the top and bottom 10%, along with their occurrences, as illustrated in Figure 3. In Figure 3a, nine cities—Huangshan, Fangchenggang, Sanya, Zhangjiajie, Bazhong, Zhoushan, Haikou, Ziyang, and Qingyang—exhibit higher industrial green TFPs, with occurrences in the top 10% exceeding 14 times. Conversely, Figure 3b highlights another nine cities—Jincheng, Datong, Taiyuan, Suzhou, Tangshan, Quanzhou, Jiaxing, Chengdu, and Harbin—that have lower industrial green TFPs, with their frequencies in the lowest 10% surpassing 12 occurrences.

3.6. Sensitivity Analysis of the Calculated Results

In Figure 3a,b, the calculated industrial green TFPs for cities such as Qingyang, Chengdu, Suzhou, and Quanzhou deviate somewhat from initial expectations. Typically, Suzhou, Quanzhou, and Chengdu are relatively affluent regions, which would suggest higher industrial green TFPs, whereas Qingyang, being less developed, might be expected to show lower TFPs. To further investigate the reasons behind the variations in urban industrial green TFPs—whether notably high or low—this paper employs Equations (11) to (14) to conduct a sensitivity analysis of the results.

$$Pro{b}_{max,yy}=\frac{{\displaystyle \sum _T{\displaystyle \sum _{Num\_Max}I\left\{exp\left(Y{Y}_{num\_max}\right)>prctile\left[exp\left(yy\right),p\right]\right\}}}}{T\times Num\_Max}$$

(11)

$$Pro{b}_{max,z{z}_{\phi }}=\frac{{\displaystyle \sum _T{\displaystyle \sum _{Num\_Max}I\left\{exp\left(Z{Z}_{\phi ,num\_max}\right)<prctile\left[exp\left(z{z}_{\phi }\right),100-p\right]\right\}}}}{T\times Num\_Max}$$

(12)

$$Pro{b}_{min,yy}=\frac{{\displaystyle \sum _T{\displaystyle \sum _{Num\_Min}I\left\{exp\left(Y{Y}_{num\_\min }\right)<prctile\left[exp\left(yy\right),q\right]\right\}}}}{T\times Num\_Min}$$

(13)

$$Pro{b}_{min,z{z}_{\phi }}=\frac{{\displaystyle \sum _T{\displaystyle \sum _{Num\_Min}I\left\{exp\left(Z{Z}_{\phi ,num\_\min }\right)>prctile\left[exp\left(z{z}_{\phi }\right),100-q\right]\right\}}}}{T\times Num\_Min}$$

(14)

In Equations from (11) to (14), $num\_max$ represents the assigned numbers for the nine cities depicted in Figure 3a, which exhibit relatively high industrial green TFPs. Conversely, $num\_min$ denotes the assigned numbers for the nine cities shown in Figure 3b that display relatively poor industrial green TFP outcomes; $Num\_Max$ and $Num\_Min$ are the total number of cities that have relatively good or poor calculated results, $Num\_Max=9$, $Num\_Min=9$; $T$ represents the total number of years covered by the research period, $T=17$; $\exp \left(·\right)$ is the power function based on natural number, while $prctile\left[·\right]$ is the percentile function, where $p$ and $q$ separately indicate the $p^{th}$ and $q^{th}$ percentiles; $I\left\{·\right\}$ is an indicator function that takes the value of 1 when the expression within the braces is true; $YY$ indicates the panel data of the desired industrial outputs per capita with the spatial units in the rows and the periods in the columns, and $yy$ is the stacked sequence of $YY$ with the spatial units pooled firstly and then the periods; $Z{Z}_{\phi }$ indicates the panel data of both the industrial inputs per capita and the share of the undesired outputs in total outputs, $Z{Z}_{\phi }=\left\{kk,ee,y_{WW},y_{s{o}_2},y_{s\&d}\right\}$, $z{z}_{\phi }$ indicates the corresponding stacked sequences. Furthermore, $Pro{b}_{max,yy}$, $Pro{b}_{max,z{z}_{\phi }}$, $Pro{b}_{min,yy}$, and $Pro{b}_{min,z{z}_{\phi }}$ are probabilities defined as in equations from (11) to (14), wherein the former two are for the cities whose calculated industrial green TFPs are good, while the latter two are for the cities whose calculated results are poor. The thresholds of $Pro{b}_{max,z{z}_{\phi }}$ and $Pro{b}_{min,z{z}_{\phi }}$ are defined separately as the $100-p$ percentile and the $100-q$ percentile of the corresponding data. Figure 4 presents the results of the sensitivity analysis.

In Figure 4a, the horizontal axis represents the value of $p$, $p\in \left[50,100\right)$, while the vertical axis indicates both $Pro{b}_{max,yy}$ and $Pro{b}_{max,z{z}_{\phi }}$; in Figure 4b, the horizontal axis is the value of $q$, $q\in \left(0,50\right]$, while the vertical axis indicates both $Pro{b}_{min,yy}$ and $Pro{b}_{min,z{z}_{\phi }}$. Figure 4a illustrates that as the quantile point decreases from 100 to 50, the probability that per capita outputs in the nine cities with superior industrial green TFP results exceed their critical values increases. Additionally, the likelihood of per capita inputs and the ratio of undesired to desired outputs falling below their critical values also rises. This indicates that with higher per capita desired outputs, lower per capita inputs, and a smaller proportion of undesired outputs relative to total outputs, the closer the critical value approaches to the median and the higher the probability of achieving better industrial green TFP outcomes. The primary factors contributing to excellent industrial green TFPs include high desired outputs per capita, low inputs per capita, and a low share of undesired outputs in total outputs. Conversely, Figure 4b shows that as the quantile point decreases from 50 to 1, the probability of the per capita outputs in the nine cities with poorer TFP results being below the critical value increases, as does the likelihood of per capita inputs and the share of undesired outputs relative to total outputs exceeding their critical values. This suggests that the principal reasons for poor industrial green TFPs include lower desired outputs per capita, higher inputs per capita, and a higher proportion of undesired outputs in total outputs.

Figure 5 elucidates the varied reasons behind the relatively high or low industrial green TFP results for individual cities. In Figure 5a–c, the vertical axis measures the frequency with which the desired outputs per capita fall below their critical values, or the inputs per capita and the shares of undesired outputs in total outputs exceed their respective critical values over the period from 2003 to 2019. The critical values on the horizontal axes are set at specific percentiles including 10, 20, 30, 40, and 50. Conversely, Figure 5d displays a vertical axis that tracks the frequency of instances where the desired outputs per capita exceed their critical value, or where the inputs per capita and the shares of undesired outputs in total outputs are below their corresponding critical values. The critical values in this figure are established at higher percentiles, specifically 55, 65, 75, 85, and 95, reflecting different thresholds for assessing performance.

In Figure 5d, the industrial sector in Qingyang shows lower shares of undesired outputs—including industrial wastewater, sulfur dioxide emissions, and smoke (dust) emissions—which contribute to its relatively high calculated green TFPs. Conversely, in Figure 5a, Suzhou’s industrial sector is characterized by lower desired outputs per capita combined with higher shares of the three undesired outputs, leading to relatively poor green TFP results. Similarly, in Figure 5b, the industrial sector in Quanzhou faces challenges akin to those in Suzhou, resulting in poor green TFPs. Additionally, Quanzhou’s relatively higher per capita industrial capital and energy inputs further contribute to its suboptimal performance. In Figure 5c, Chengdu exhibits poor green TFPs for reasons similar to those observed in Quanzhou. However, the share of smoke (dust) emissions in Chengdu’s total outputs does not significantly impact its green TFP assessment, indicating that other factors primarily drive its poor performance.

4. Assessing the Accuracy of STE-SRM: Comparative Analysis with DEA-SBM and Bayesian SFA

Data Envelopment Analysis (DEA) and Stochastic Frontier Analysis (SFA) represent two additional methodological approaches for calculating green Total Factor Productivities (TFPs). This paper recalculates the industrial green TFPs for over 280 cities from 2003 to 2019 using the DEA-SBM from the DEA family, and subsequently through the Bayesian SFA from the SFA family. It then assesses the accuracy of the STE-SRM in calculating green TFPs by comparing these results.

4.1. Reassessment of China’s Urban Industrial Green TFPs Using DEA-SBM

As an important method within the DEA family, the DEA-SBM is widely used in the calculation of green Total Factor Productivities (TFPs). This model thoroughly considers the slackness of input and output variables and evaluates changes in inputs and outputs from both non-radial and non-angular perspectives simultaneously. The main features of the DEA-SBM model include the following: firstly, the slackness of inputs and outputs is incorporated into the objective function; secondly, changes in inputs and outputs are examined simultaneously; thirdly, it allows for the examination of different proportional changes in inputs and outputs; fourthly, undesirable outputs can be included in the efficiency calculation.

We consider each city as an industrial production decision-making unit. We set $i=1,2,\cdots ,N$ as the number of the cities, and define $x_{k,i}$, $y_{r,i}$, and $y_{\pi ,i}$ as the $k^{th}$ input, the $r^{th}$ desired output, and the ${\pi }^{th}$ undesired output of the city $i$ respectively, wherein $k=1,2,\cdots ,m_0$, $r=1,2,\cdots ,m_1$, and $\pi =1,2,\cdots ,m_2$, then the production possibility sets of the industrial sectors in the city $i$ can be defined as in the following Equation (15).

$$P_i=\left\{\left.\left(x_{k,i},y_{r,i},y_{\pi ,i}\right)\right|x_{k,i}\ge {\displaystyle \sum _l{\lambda }_l{x}_{k,l}},y_{r,i}\le {\displaystyle \sum _l{\lambda }_l{y}_{r,l}},y_{\pi ,i}\ge {\displaystyle \sum _l{\lambda }_l{y}_{\pi ,l}},{\lambda }_l\ge 0\right\}$$

(15)

As detailed in Section 3.2, let us assume that within the industrial production decision-making framework, the inputs consist of industrial capital, labor, and energy. The outputs are categorized into one desired and three undesired types. The desired output is the gross value of industrial production, whereas the undesired outputs encompass industrial wastewater emissions, sulfur dioxide emissions, and smoke (dust) emissions. Thus, in Equation (15), $m_0=3$, $m_1=1$, and $m_2=3$. In Equation (15), $l$ also represents the number of cities, $l=1,2,\cdots ,N$; $P_i$ represents the production possibility sets of the industrial production of the city $i$; ${\lambda }_l$ represents the weight of the city $l$. According to the DEA-SBM, we define $s_k^{-}$, $s_r^{+}$, and $s_{\pi }^{-}$ as the slacks of the $k^{th}$ input, the $r^{th}$ desired output, and the ${\pi }^{th}$ undesired output, respectively; the green TFPs can be calculated by the following two Equations (16) and (17).

$$GTF{P}_{dea,i}=Min\left\{\frac{1-\frac{1}{m_0}{\displaystyle \sum _{k=1}^{m_0}\left(\frac{s_k^{\_}}{x_{k,i}}\right)}}{1+\frac{1}{m_1+m_2}\left[{\displaystyle \sum _{r=1}^{m_1}\left(\frac{s_r^{\_}}{y_{r,i}}\right)}+{\displaystyle \sum _{\pi =1}^{m_2}\left(\frac{s_{\pi }^{+}}{y_{\pi ,i}}\right)}\right]}\right\}$$

(16)

$$\begin{array}{l}s.t.\\ \left\{\begin{array}{l}{x}_{k,i}={\displaystyle \sum _l{\lambda }_l{x}_{k,l}}+s_k^{-}\\ y_{r,i}={\displaystyle \sum _l{\lambda }_l{y}_{r,l}}-s_r^{+}\\ y_{\pi ,i}={\displaystyle \sum _l{\lambda }_l{y}_{\pi ,l}}+s_{\pi }^{-}\\ \lambda \ge 0,s_k^{-}\ge 0,s_r^{+}\ge 0,s_{\pi }^{-}\ge 0\end{array}\right.\end{array}$$

(17)

In Equation (16), $GTF{P}_{dea,i}$ is the industrial green TFP of the city $i$ calculated by DEA-SBM, and it is the decreasing function of $s_k^{-}$, $s_r^{+}$, and $s_{\pi }^{-}$. If the values of $s_k^{-}$, $s_r^{+}$, and $s_{\pi }^{-}$ are zeros, the value of $GTF{P}_{dea,i}$ is biggest, $GTF{P}_{dea,i}=1$. If the values of $s_k^{-}$, $s_r^{+}$, and $s_{\pi }^{-}$ are respectively equal to $x_{k,i}$, $y_{r,i}$, and $y_{\pi ,i}$, the value of $GTF{P}_{dea,i}$ will be lowest, $GTF{P}_{dea,i}=0$. Thus, $GTF{P}_{dea,i}\in \left[0,1\right]$. Based on Equations (16) and (17), and utilizing MATLAB R2023a along with custom codes developed by the author, the industrial green Total Factor Productivities (TFPs) of 280 Chinese cities from 2003 to 2019 have been recalculated, as illustrated in Figure 6b.

4.2. Revised Calculation of China’s Urban Industrial Green TFPs Using Bayesian SFA

4.2.1. Fundamental Logic of Bayesian SFA in Calculating Green TFPs

Bayesian SFA is a prominent method within the SFA family, representing an integration of stochastic frontier analysis with Bayesian methods. In the assessment of green production efficiency, Bayesian SFA primarily focuses on determining the posterior distribution of unknown parameters in the stochastic frontier production function model. This is achieved by leveraging their prior distributions, in conjunction with the total probability formula and sample data. The green production efficiency is then calculated using these estimated parameters. We define $Y$ as the outputs and $Z_{\upsilon }$ as the ${\upsilon }^{th}$ input, $\upsilon =1,2,\cdots ,\Upsilon$, and set $y^{\ast }=Ln\left(Y\right)$, $z_v^{\ast }=Ln\left(Z_{\upsilon }\right)$, and $z^{\ast }=\left\{z_v^{\ast }\right\}$; the stochastic frontier production function model can be preset as the following Equation (18).

$$y^{\ast }=g\left(z^{\ast };\gamma \right)+\xi -\varsigma$$

(18)

In Equation (18), $\gamma$ represents the exogenous parameters of the inputs and their various combinations; $\xi$ is a disturbance term that follows a normal distribution with a zero mean and a specific variance, $\xi ~N\left(0,{\sigma }_{\xi }^2{I}_{NT}\right)$; the indices $i$ and $t$ represent the analyzed regions and periods, respectively, $i=1,2,\cdots ,N$, $t=1,2,\cdots ,T$; $\varsigma$ is the non-efficiency term and subject to the semi-normal distribution, $\varsigma ~N^{+}\left(0,{\sigma }_{\varsigma }^2{I}_{NT}\right)$; $\varsigma$ and $\xi$ are independent from each other. In Bayesian SFA, the green TFPs are usually defined by $GTF{P}_{sfa,it}=\exp \left(-{\varsigma }_{it}\right)$, and their values usually range between 0 and 1. In Equation (18), $g\left(z^{\ast };\gamma \right)$ can vary depending on the different input factors involved. Drawing on the analyses by Makiela (2014) and Makiela and Ouattara (2018) [36,37], the general formula for $g\left(z^{\ast };\gamma \right)$ is established as Equation (19).

$$g\left(z^{\ast };\gamma \right)={\displaystyle \sum _{\kappa =0}^{\Lambda }{t}^{\kappa }{\gamma }_{0,\kappa }}+{\displaystyle \sum _{\kappa =0}^{\Lambda }\left\{{\displaystyle \sum _{\upsilon =1}^{\Upsilon }\left(t^{\kappa }{z}_{\upsilon }^{\ast }\right)}{\gamma }_{{\upsilon }_1}+{\displaystyle \sum _{\upsilon =1}^{\Upsilon }\left[t^{\kappa }{\left(z_{\upsilon }^{\ast }\right)}^2\right]}{\gamma }_{{\upsilon }_2}+{\displaystyle \sum _{\tau =1}^{\Upsilon }{\displaystyle \sum _{\upsilon =1}^{\Upsilon }\left(t^{\kappa }{z}_{\tau }^{\ast }{z}_{\upsilon }^{\ast }\right){\gamma }_{{\upsilon }_3}}}\right\}}$$

(19)

In Equation (19), $\upsilon$ denotes the subscript for the inputs, with $\upsilon =1,2,\cdots ,\Upsilon$; the symbols ${\upsilon }_1$, ${\upsilon }_2$, and ${\upsilon }_3$ are the subscripts for the exogenous parameters, where ${\upsilon }_1=\upsilon +\kappa C_{\Upsilon }^2+2\kappa \Upsilon$, ${\upsilon }_2=\upsilon +\kappa C_{\Upsilon }^2+\left(2\kappa +1\right)\Upsilon$, and ${\upsilon }_3=\psi +\kappa C_{\Upsilon }^2+\left(2\kappa +2\right)\Upsilon$. Wherein, $\kappa$ indicates the exponential value of the period variable $t$, $\kappa =0,1,\cdots ,\Lambda$, and $\Lambda$ can be 0, 1, or 2; $\psi$ represents the number of exogenous parameters associated with the new variables created by the cross-multiplication of two inputs, $\psi =1,2,\cdots ,C_{\Upsilon }^2$; $C_{\Upsilon }^2$ is the total number of the new variables formed, $C_{\Upsilon }^2=\frac{\Upsilon !}{\left(\Upsilon -2\right)!2!}$.

We define $g\left(z^{\ast };\gamma \right)=z^{\ast }\gamma$, and $\gamma$ subjects to a normal distribution, $\gamma ~N\left(b,C^{-1}\right)$. Also, we respectively define ${\sigma }_{\xi }^2$ and ${\sigma }_{\varsigma }^2$ subject to the inverse gamma distribution as ${\sigma }_{\xi }^2~IG\left(\eta n_0,\eta {\alpha }_0\right)$, ${\sigma }_{\varsigma }^2~IG\left(5,10L{n}^2\left(GTF{P}_{sfa}^0\right)\right)$, where $\eta$ represents an empirical constant, $\eta =0.5$, and $GTF{P}_{sfa}^0$ is the prior mean of the green TFPs. Moreover, we define $f_N\left(\left.·\right|△,\square \right)$ as representing the probability density function of the normal distribution with the mean of $△$ and the variance of $\square$, and define $f_G\left(\left.·\right|△,\square \right)$ as representing the probability density function of the gamma distribution with the mean of $△/\square$ and the variance of ${\left(△/\square \right)}^2$; then, Equation (18) can be changed to its Bayesian form, as shown in Equation (20).

$$f_N\left(\left.\gamma \right|b,C^{-1}\right)·f_G\left(\left.{\sigma }_{\xi }^{-2}\right|\eta n_0,\eta {\alpha }_0\right)·f_G\left(\left.{\sigma }_{\varsigma }^{-2}\right|5,10L{n}^2\left(GTF{P}_{sfa}^0\right)\right)·{\displaystyle \prod _{t=1}^T{\displaystyle \prod _{i=1}^N{f}_N\left(\left.y_{it}^{\ast }\right|z_{it}^{\ast }\gamma -{\varsigma }_{it},{\sigma }_{\xi }^2\right)f_G\left(\left.{\varsigma }_{it}\right|0,{\sigma }_{\varsigma }^2\right)}}$$

(20)

According to the Bayesian analysis, the conditional posterior distribution of the unknown parameters of $\gamma$, ${\sigma }_{\xi }^2$, ${\sigma }_{\varsigma }^2$, and ${\varsigma }_{it}$ can be shown as Equations (21) to (24).

$$p\left(\left.\gamma \right|y^{\ast },z^{\ast },{\varsigma }_{it},{\sigma }_{\xi }^{-2},{\sigma }_{\varsigma }^{-2}\right)\propto f_N^J\left(\left.\gamma \right|C_{\ast }^{-1}\left[Cb+{\sigma }_{\xi }^{-2}{z}^{\ast \prime }\left(y^{\ast }+{\varsigma }_{it}\right)\right],C_{\ast }^{-1}\right)$$

(21)

$$p\left(\left.{\sigma }_{\xi }^{-2}\right|y^{\ast },z^{\ast },{\varsigma }_{it},{\sigma }_{\varsigma }^{-2},\gamma \right)\propto f_G\left(\left.{\sigma }_{\xi }^{-2}\right|\frac{n_0+NT}{2},\frac{a_0+{\xi }^{\prime }\xi }{2}\right)$$

(22)

$$p\left(\left.{\sigma }_{\varsigma }^{-2}\right|y^{\ast },z^{\ast },{\varsigma }_{it},{\sigma }_{\xi }^{-2},\gamma \right)\propto f_G\left(\left.{\sigma }_{\varsigma }^{-2}\right|\frac{NT}{2}+5,\frac{1}{2}{\displaystyle \sum _{t=1}^T{\displaystyle \sum _{i=1}^N{\varsigma }_{it}}}+10L{n}^2\left(GTF{P}_{sfa}^0\right)\right)$$

(23)

$$p\left(\left.{\varsigma }_{it}\right|y^{\ast },z^{\ast },{\sigma }_{\xi }^{-2},{\sigma }_{\varsigma }^{-2},\gamma \right)\propto f_N^{NT}\left(\left.{\varsigma }_{it}\right|\frac{{\sigma }_{\varsigma }\left(z^{\ast }\gamma -y^{\ast }\right)}{{\sigma }_{\varsigma }^2+{\sigma }_{\xi }^2},\frac{{\sigma }_{\varsigma }^2{\sigma }_{\xi }^2}{{\sigma }_{\varsigma }^2+{\sigma }_{\xi }^2}\right)I\left({\varsigma }_{it}\in R_{+}^{NT}\right)$$

(24)

In Equations (21) and (22), $C_{\ast }^{-1}={\left(C+{\sigma }_{\xi }^{-2}{z}^{\ast \prime }{z}^{\ast }\right)}^{-1}$, $\xi =y^{\ast }+{\varsigma }_{it}-z^{\ast }\gamma$, $I\left({\varsigma }_{it}\in R_{+}^{NT}\right)$ is an indicative function indicator that $\varsigma$ obeys the semi-normal distribution. Based on Equations (20) to (24) and employing techniques such as Gibbs sampling and Metropolis-Hastings sampling [38], the unknown parameters can be effectively estimated and inferred through posterior analysis.

4.2.2. China’s Urban Industrial Green TFPs Re-calculated by the Bayesian SFA

The general term formula in Equation (19) is somewhat complex and may lead to multicollinearity issues when modeling China’s urban industrial production process. Therefore, in the recalibration of industrial green TFPs for China’s 280 cities from 2003 to 2019 using Bayesian SFA, Equation (19) has been constrained to Equation (25) to mitigate these challenges.

$$g^{″}\left(z^{\ast ″};\gamma \right)={\gamma }_{0,0}+{\displaystyle \sum _{\upsilon =1}^6{z}_{\upsilon }^{\ast ″}}{\gamma }_{\upsilon }+{\displaystyle \sum _{\upsilon =1}^6{\left(z_{\upsilon }^{\ast ″}\right)}^2{\gamma }_{\upsilon +6}}+{\displaystyle \sum _{\tau =1}^6{\displaystyle \sum _{\upsilon =1}^6\left(z_{\tau }^{\ast ″}{z}_{\upsilon }^{\ast ″}\right){\gamma }_{\psi +12}}}+{\gamma }_{0,1}t+{\displaystyle \sum _{\upsilon =1}^6\left(t{z}_{\upsilon }^{\ast ″}\right)}{\gamma }_{\upsilon +27}$$

(25)

In Equation (25), six inputs are considered, $z^{\ast ″}=\left\{k^{\ast },l^{\ast },e{n}^{\ast },y_{ww}^{\ast },y_{s{o}_2}^{\ast },y_{s\&d}^{\ast }\right\}$; these inputs include three real inputs of capital, labor, and energy, and three undesired outputs. Wherein, $k^{\ast }=Ln\left(K\right)$, $l^{\ast }=Ln\left(L\right)$, $e{n}^{\ast }=Ln\left(En\right)$, $y_{ww}^{\ast }=-Ln\left(Y_{WW}\right)$, $y_{s{o}_2}^{\ast }=-Ln\left(Y_{S{O}_2}\right)$, $y_{s\&d}^{\ast }=-Ln\left(Y_{S\&D}\right)$, and the definitions of $K$, $L$, $En$, $Y_{WW}$, $Y_{S{O}_2}$, and $Y_{S\&D}$ are the same as before. At this time, $\Upsilon =6$. Meanwhile, if $\Lambda$ is preset to be 1, then $\kappa =0,1$, $C_{\Upsilon }^2=15$, and $\psi =1,2,\cdots ,15$. It is important to note that the revised formula in Equation (20) omits the interaction term between the period variable $t$ and the quadratic forms of the inputs. Additionally, it does not include terms representing the product of $t$ and the cross-multiplication of inputs paired together.

From Equation (25), five distinct degradation models can be derived, each based on different underlying assumptions. There will be a basic Cobb–Douglas model (C-D) if ${\gamma }_{\upsilon +6}=0$, ${\gamma }_{\psi +12}=0$, ${\gamma }_{0,1}=0$, and ${\gamma }_{\upsilon +27}=0$, a Cobb–Douglas model with the period variable (CD-T) if ${\gamma }_{\upsilon +6}=0$, ${\gamma }_{\psi +12}=0$, and ${\gamma }_{\upsilon +27}=0$, and a Cobb–Douglas model with the parameters varied with the time (CD-Linear-T) if ${\gamma }_{\upsilon +6}=0$ and ${\gamma }_{\psi +12}=0$. Meanwhile, there will be a basic transcendental logarithmic model (Trans-Log) if ${\gamma }_{0,1}=0$ and ${\gamma }_{\upsilon +27}=0$, and a Trans-Log model with the period variable (Trans-Log-T) if ${\gamma }_{\upsilon +27}=0$. Based on theoretical process similar to the equations from (20) to (24), these five degradation models can be estimated, as shown in

Table 4. Bayesian estimators of different stochastic frontier production functions of China’s urban industrial sectors.

	C-D	CD-T	CD-Linear-T	Trans-Log	Trans-Log-T
$Const.$	0.6462 (57.82) ***	0.5184 (40.24) ***	−0.0383 (−16.75) ***	−0.1203 (−10.96) ***	−0.0997 (−9.29) ***
$Ln\left(K\right)$	0.3437 (26.83) ***	0.4055 (32.11) ***	0.0133 (5.48) ***	0.0408 (2.09) **	0.0221 (1.17)
$Ln\left(L\right)$	0.1136 (10.57) ***	0.0676 (6.34) ***	0.0014 (0.68)	0.0431 (3.15) ***	0.0139 (1.03)
$Ln\left(En\right)$	−0.1846 (−19.07) ***	−0.2218 (−23.37) ***	−0.0076 (−4.08) ***	−0.0618 (−4.58) ***	−0.0606 (−4.66) ***
$Ln\left(y_{ww}\right)$	−0.0314 (−3.34) ***	0.0246 (2.62) ***	−0.0085 (−4.66) ***	−0.0058 (−0.44)	0.0147 (1.14)
$Ln\left(y_{S{O}_2}\right)$	−0.0999 (−12.19) ***	−0.0824 (−10.43) ***	0.0076 (4.74) ***	0.0377 (3.35) ***	0.0220 (2.01) **
$Ln\left(y_{S\&D}\right)$	−0.4806 (−6.25) ***	−0.0417 (−18.23) ***	−0.1344 (−9.89) ***	−0.9893 (−9.28) ***	−0.8297 (−8.06) ***
${\left[Ln\left(K\right)\right]}^2$				−0.0457 (−3.62) ***	−0.0271 (−2.17) **
${\left[Ln\left(L\right)\right]}^2$				0.0776 (4.86) ***	0.0722 (4.64) ***
${\left[Ln\left(En\right)\right]}^2$				0.0145 (0.87)	0.0119 (0.72)
${\left[Ln\left(y_{ww}\right)\right]}^2$				0.0202 (1.35)	0.0277 (1.92) **
${\left[Ln\left(y_{s_2}\right)\right]}^2$				0.0435 (3.45) ***	0.0401 (3.27) ***
${\left[Ln\left(y_{s\&d}\right)\right]}^2$				0.1826 (1.46)	0.2570 (2.10) **
$Ln\left(K\right)\times Ln\left(L\right)$				−0.0406 (−6.17) ***	−0.0277 (−4.26) ***
$Ln\left(K\right)\times Ln\left(En\right)$				−0.0215 (−1.62)	−0.0221 (−1.73) *
$Ln\left(K\right)\times Ln\left(y_{ww}\right)$				−0.0479 (−4.18) ***	−0.0281 (−2.51) **
$Ln\left(K\right)\times Ln\left(y_{S{O}_2}\right)$				0.0421 (3.72) ***	0.0345 (3.11) ***
$Ln\left(K\right)\times Ln\left(y_{S\&D}\right)$				0.0543 (0.59)	0.0723 (0.82)
$Ln\left(L\right)\times Ln\left(En\right)$				0.0536 (6.78) ***	−0.0534 (−6.82) ***
$Ln\left(L\right)\times Ln\left(y_{ww}\right)$				−0.0881 (−7.71) ***	−0.0843 (−7.51) ***
$Ln\left(L\right)\times Ln\left(y_{S{O}_2}\right)$				0.0116 (1.18)	0.0150 (1.57)
$Ln\left(L\right)\times Ln\left(y_{S\&D}\right)$				0.6164 (5.78) ***	0.6585 (6.23) ***
$Ln\left(En\right)\times Ln\left(y_{ww}\right)$				0.0300 (5.22) ***	0.0257 (4.61) ***
$Ln\left(En\right)\times Ln\left(y_{S{O}_2}\right)$				−0.0244 (−2.92) ***	−0.0190 (−2.32) ***
$Ln\left(En\right)\times Ln\left(y_{S\&D}\right)$				−0.4508 (−5.54) ***	−0.5236 (−6.59) ***
$Ln\left(y_{ww}\right)\times Ln\left(y_{S{O}_2}\right)$				−0.0109 (−3.02) ***	−0.0059 (−1.67) *
$Ln\left(y_{ww}\right)\times Ln\left(y_{S\&D}\right)$				0.2608 (3.66) ***	0.1262 (1.81) *
$Ln\left(y_{S{O}_2}\right)\times Ln\left(y_{S\&D}\right)$				−1.1460 (−2.73) ***	0.0396 (16.15) ***
$T$		0.3842 (5.33) ***	0.8187 (37.88) ***		−1.2758 (−3.09) ***
$T\times Ln\left(K\right)$			0.3084 (12.09) ***
$T\times Ln\left(L\right)$			0.0655 (3.41) ***
$T\times Ln\left(En\right)$			−0.1572 (−8.62) ***
$T\times Ln\left(y_{ww}\right)$			0.0553 (3.02) ***
$T\times Ln\left(y_{S{O}_2}\right)$			−0.1244 (−8.58) ***
$T\times Ln\left(y_{S\&D}\right)$			0.2432 (1.83) *
${\sigma }_{\varsigma }^2$	0.3458 (12.94) ***	0.3176 (13.59) ***	0.2923 (13.94) ***	0.3070 (12.38) ***	0.2876 (12.29) ***
${\sigma }_{\xi }$	0.3562 (29.69) ***	0.3475 (32.51) ***	0.3307 (33.35) ***	0.3460 (29.81) ***	0.3387 (30.54) ***
MDD	−1050.55	−975.72	−862.52	−1024.41	−967.64
SSE	322.88	311.53	280.00	310.59	300.14

Note: the above outputs were collected based on MATLAB R2023a. ( ) represents T-statistic, ***, ** and * mean having passed the hypothesis test with a significance level of 1%, 5% and 10%. MDD indicates marginal data density based on Harmonic Mean Estimator, while SSE indicates the sum of squared errors in the Batesian SFA model.

. It is worth noting that, in the Bayesian sampling process, the prior mean of China’s urban industrial green TFPs is preset to be 0.75, $GTF{P}_{sfa}^0=0.75$; $n_0$ and ${\alpha }_0$ are both preset to be 10⁻⁶, $n_0={\alpha }_0={10}^{-6}$; the total number of samples taken is 50,000. The parameters are estimated and inferred based on the mean of these parameters over the last 40,000 samples, after discarding the initial 10,000 samples.

The estimation results presented in Table 4 reveal that several parameters for the explanatory variables in the CD-Linear-T, Trans-Log, and Trans-Log-T models are not statistically significant, suggesting that these models may not be optimal for China’s urban industrial sectors. Additionally, in the CD-T model, the elasticity of wastewater to outputs is positively correlated, which contradicts expected outcomes, thereby disqualifying CD-T as the optimal model as well. Consequently, the Cobb–Douglas (C-D) model remains for use in the Bayesian SFA to calculate the industrial green TFPs. In practice, the C-D model is considered optimal because it exhibits statistically significant parameters and robust statistical characteristics. From the estimated technical efficiency items of the C-D model in Table 4, combined with the formula of $GTF{P}_{sfa,it}=\exp \left(-{\varsigma }_{it}\right)$, China’s urban industrial green TFPs can be re-calculated by the Bayesian SFA, as shown in Figure 6c.

4.2.3. Comparative Analysis of Three Methods and Accuracy Assessment of the STE-SRM

To evaluate the accuracy of the STE-SRM, this section compares the calculated results obtained using the three different methods. As illustrated in Figure 6, the industrial green TFPs for China’s 280 cities from 2003 to 2019 are depicted in subgraphs 6a, 6b, and 6c. Additionally, lognormal distribution fitting of these results is shown in subgraph 6d. From subgraph 6a, the range of urban industrial green TFPs calculated by the STE-SRM spans from 0 to 40, whereas in subgraphs 6b and 6c, the TFPs calculated by the DEA-SBM and Bayesian SFA range from 0 to 1. Subgraph 6d shows that the cumulative density functions of all three types of industrial green TFPs follow power function trends, where the exponents for the TFPs calculated by the DEA-SBM and Bayesian SFA exceed one, while the exponent for the TFPs calculated by the STE-SRM is less than one. Moreover, the power function curve for the results from Bayesian SFA is steeper than that from the DEA-SBM.

Additionally, the main statistical characteristics of the three types of industrial green TFPs are detailed in

Table 5. Key descriptive statistics of industrial green TFPs derived from various methods.

	STE-SRM	Bayesian SFA	DEA-SBM
Mean	0.9872	0.6628	0.4640
Median	0.6260	0.6933	0.4308
Maximum	38.2525	0.9407	1.0000
Minimum	0.0663	0.1533	0.0021
Std. Dev.	1.4340	0.1404	0.2250
Coefficient of Variation	0.6884	4.7208	2.0622
75th Quartile	1.0893	0.7698	0.5769
25th Quartile	0.3729	0.5768	0.3080
Quartile Deviation	0.7163	0.1931	0.2689

Note: the results were obtained based on MATLAB R2023a and EViews 11.0.

. According to the table, the STE-SRM calculated values have the highest mean, maximum, and standard deviation, with the median and minimum values falling in between. Furthermore, the quartile deviation (the difference between the 75th and 25th quartiles) of the STE-SRM results is the largest, and its coefficient of variation (mean divided by the standard deviation) is the smallest, indicating a wider range but higher concentration of the TFPs calculated by the STE-SRM. Conversely, the Bayesian SFA results show the largest coefficient of variation and the smallest quartile deviation, suggesting a narrower range but greater dispersion. The DEA-SBM’s results fall between these extremes in both coefficient of variation and quartile deviation.

In summary, although the theoretical frameworks for calculating industrial green TFPs using STE-SRM, DEA-SBM, and Bayesian SFA differ and lack direct comparability, the comparative analysis of the results suggests that the STE-SRM is suitable for researchers preferring a broader range with concentrated outcomes, while the Bayesian SFA is better for those seeking a narrower range with more dispersed results.

5. Advanced Analysis of China’s Urban Industrial Green TFPs Calculated by STE-SRM: Examining Spatial Heterogeneity and Convergence

5.1. Spatial Heterogeneity Analyzed Using Dagum’s Gini Coefficient

Dagum’s Gini coefficient is a significant tool for analyzing spatial heterogeneity. Its fundamental principle involves decomposing the total coefficient into three distinct parts: the contribution from intra-group differences, the contribution from inter-group differences, and the contribution from the intensity of trans-variation differences [39]. This section explores the spatial heterogeneity of China’s urban industrial green TFPs using Dagum’s Gini coefficient, as detailed in Equation (26). Here, $G_{GTFP}$ represents the total Dagum’s Gini coefficient, while $G_{GTFP,w}$, $G_{GTFP,nb}$, and $G_{GTFP,t}$ denote the three decomposed components, respectively.

$$G_{GTFP}=\underset{G_{GTFP,w}}{\underbrace{{\displaystyle \sum _{h=1}^k{G}_{hh}{P}_h{S}_h}}}+\underset{G_{GTFP,nb}}{\underbrace{{\displaystyle \sum _{h=1}^k{\displaystyle \sum _{h^{\prime }\ne h}{G}_{h{h}^{\prime }}\left(P_h{S}_{h^{\prime }}+P_{h^{\prime }}{S}_h\right)}}{D}_{h{h}^{\prime }}}}+\underset{G_{GTFP,t}}{\underbrace{{\displaystyle \sum _{h=1}^k{\displaystyle \sum _{h^{\prime }\ne h}{G}_{h{h}^{\prime }}\left(P_h{S}_{h^{\prime }}+P_{h^{\prime }}{S}_h\right)}}\left(1-D_{h{h}^{\prime }}\right)}}$$

(26)

In Equation (26), $G_{hh}$ and $G_{h{h}^{\prime }}$ respectively indicate the intra-group Gini coefficient and the inter-group Gini coefficient, $G_{hh}=\frac{{\displaystyle \sum _{i=1}^{n_h}{\displaystyle \sum _{l=1}^{n_h}\left|\widehat{GTF{P}_{hi}}-\widehat{GTF{P}_{hl}}\right|}}}{2n_h^2\overline{\widehat{GTF{P}_h}}}$, $G_{h{h}^{\prime }}=\frac{{\displaystyle \sum _{i=1}^{n_h}{\displaystyle \sum _{r=1}^{n_{h^{\prime }}}\left|\widehat{GTF{P}_{hi}}-\widehat{GTF{P}_{h^{\prime }r}}\right|}}}{n_h{n}_{h^{\prime }}\left(\overline{\widehat{GTF{P}_h}}+\overline{\widehat{GTF{P}_{h^{\prime }}}}\right)}$; $h$ and $h^{\prime }$ respectively represent the spatial groups, where in this paper, the spatial groups are classified by eastern, central, western, and northeast China, and then $h=1,2,3,4$, $h^{\prime }=1,2,3,4$; $P_h$ indicates the ratios of the number of the cities belong to the spatial group $h$ to the total number of the cities, $P_h=N_h/N$; $S_h$ and $S_{h^{\prime }}$ respectively indicate the ratios of the sum of the industrial green TFPs of the cities belonging to the spatial group $h$ (or $h^{\prime }$) to the sum of the industrial green TFPs of all the cities, $S_h=\frac{N_h\overline{\widehat{GTF{P}_h}}}{N\overline{\widehat{GTFP}}}$, $S_{h^{\prime }}=\frac{N_{h^{\prime }}\overline{\widehat{GTF{P}_{h^{\prime }}}}}{N\widehat{GTFP}}$; $\overline{\widehat{GTF{P}_h}}$, $\overline{\widehat{GTF{P}_{h^{\prime }}}}$, and $\overline{\widehat{GTFP}}$ respectively represent the mean value of the industrial green TFPs of the spatial group $h$, the spatial group $h^{\prime }$, and all the cities; $N_h$ and $N_{h^{\prime }}$ indicate the number of the cities belonging to the spatial group $h$ and the spatial group $h^{\prime }$, respectively; $D_{h{h}^{\prime }}$ indicates the relative spatial influence between the spatial group $h$ and the spatial group $h^{\prime }$, $D_{h{h}^{\prime }}=\frac{d_{h{h}^{\prime }}-p_{h{h}^{\prime }}}{d_{h{h}^{\prime }}+p_{h{h}^{\prime }}}$; $d_{h{h}^{\prime }}$ indicates the total spatial influences of the spatial group $h$ and the spatial group $h^{\prime }$, as shown in Equation (27); $p_{h{h}^{\prime }}$ is the first moment of the trans-variation, as shown in Equation (28).

$$d_{h{h}^{\prime }}={\displaystyle \underset{0}{\overset{\infty }{\int }}d{F}_h\left(\widehat{GTF{P}_{hi}}\right){\displaystyle \underset{0}{\overset{\widehat{GTF{P}_{hi}}}{\int }}\left(\widehat{GTF{P}_{hi}}-\widehat{GTF{P}_{h^{\prime }r}}\right)d{F}_{h^{\prime }}\left(\widehat{GTF{P}_{h^{\prime }r}}\right)}}$$

(27)

$$p_{h{h}^{\prime }}={\displaystyle \underset{0}{\overset{\infty }{\int }}d{F}_{h^{\prime }}\left(\widehat{GTF{P}_{hi}}\right){\displaystyle \underset{0}{\overset{\widehat{GTFP}}{\int }}\left(\widehat{GTF{P}_{hi}}-\widehat{GTF{P}_{h^{\prime }r}}\right)d{F}_h\left(\widehat{GTF{P}_{h^{\prime }r}}\right)}}$$

(28)

Based on Equation (26) and categorizing the spatial groups by eastern, central, western, and northeast regions of China, Dagum’s Gini coefficients have been calculated and are presented in Figure 7. The division of cities includes 83 in eastern, 80 in central, 83 in western, and 34 in northeast China.

From Figure 7, the overall level of Dagum’s Gini coefficient for China’s urban industrial green TFPs ranges from approximately 0.4 to 0.5, indicating significant spatial heterogeneity. Among the four decomposed parts, the contribution from the intensity of trans-variation differences is the highest. The development trend of intra-group differences is similar to that of inter-group differences; however, the value of the intra-group differences is slightly greater than that of the inter-group differences. Regarding the trend of the overall Dagum’s Gini coefficient, the spatial heterogeneity of China’s urban industrial green TFPs increased from 2003 to 2010, followed by a general trend of narrowing after 2011. The narrowing trend decreased rapidly between 2011 and 2013, experienced some expansion from 2013 to 2016, and the expansion was subsequently controlled after 2017.

5.2. Analyzing Spatial Convergence through Sigma and Beta Convergence Methods

5.2.1. The Sigma Convergence and the Beta Convergence

To better understand the spatial attributes of China’s urban industrial green TFPs, this section employs spatial sigma and beta convergence analyses to investigate spatial convergence. In general, sigma convergence analysis assesses the coefficient of variation, as illustrated in Equation (29). A sigma convergence trend is indicated when the coefficient of variation for a specific economic indicator consistently decreases over time. In Equation (29), $Astg\_GTFP$ indicates the sigma convergence of China’s urban industrial green TFPs, $Std\left(·\right)$ and $Mean\left(·\right)$ indicate the standard deviation and the mean value, respectively, $i$ is the number of the cities, and $t$ indicates the year, while $Num\_g$ indicates the number of the cities in a special spatial group.

$$Astg\_GTF{P}_{t,Num\_g}=\underset{i\in Num\_g}{Std}\left(GTF{P}_{i,t}\right)/\underset{i\in Num\_g}{Mean}\left(GTF{P}_{i,t}\right)$$

(29)

Beta convergence is typically assessed by examining the elasticity of a specific economic indicator’s current value to its increment. The primary approach involves establishing a regression model with the logarithm of the indicator’s increment as the dependent variable and the logarithm of the current value of the indicator as the core explanatory variable. A beta convergence trend is indicated if the estimated parameter of the core explanatory variable is significantly negative. It is important to note that convergence can occur in two forms: absolute convergence and conditional convergence. Convergence is conditional if other control variables are included in the regression model; otherwise, it is absolute convergence. Moreover, if the regression model includes only the core explanatory variable and its spatiotemporal spillover terms, it also indicates absolute convergence.

Given the challenges in obtaining data on China’s industrial development across 280 cities during 2003–2019, this paper establishes a regression model, as outlined in Equations (30) and (31), to examine the spatial absolute convergence of China’s urban industrial green TFPs.

$$\begin{array}{ll}Ln\left(\frac{GTF{P}_{i,t+1}}{GTF{P}_{i,t}}\right)& ={\omega }_0\left[ST{W}_1\times Ln\left(\frac{GTF{P}_{i,t+1}}{GTF{P}_{i,t}}\right)\right]+{\omega }_{1}Ln\left(GTF{P}_{i,t}\right)\\ & +{\omega }_2\left[ST{W}_1\times Ln\left(GTF{P}_{i,t}\right)\right]+{\vartheta }_1+{\vartheta }_2+{\mu }_3\end{array}$$

(30)

$${\mu }_3={\omega }_3\left(ST{W}_1\times {\mu }_3\right)+{\epsilon }_4$$

(31)

In Equations (30) and (31), $GTFP$ still indicates the industrial green TFPs in China’s 280 cities, $t$ indicates the year, and $t+1$ indicates the lag year of $t$. $ST{W}_1$ are the spatiotemporal weight matrices constructed similarly to $STW$ previously. However, for the analyses in different spatial groups, $ST{W}_1$ should be constructed by the latitude and longitude distances among the corresponding cities belonging to the spatial group, and their temporal weight matrices are not determined by the disturbance terms in Equation (10) but determined by the different estimated residuals in the univariate linear regression models, taking $Ln\left(GTF{P}_{i,t+1}/GTF{P}_{i,t}\right)$ as the explained variable and $Ln\left(GTF{P}_{i,t}\right)$ as the explanatory variable. These univariate linear regression models also contain different sample cities when different spatial groups are considered.

In Equations (30) and (31), ${\omega }_0$ and ${\omega }_3$ are spatial correlation coefficients; ${\omega }_1$ and ${\omega }_2$ are exogenous parameters; ${\vartheta }_1$ and ${\vartheta }_2$ are parameters representing individual or period effects; ${\mu }_3$ and ${\epsilon }_4$ are disturbance terms, ${\epsilon }_4~IID\left(0,{\sigma }_4^2\right)$, and the distribution of ${\mu }_3$ is decided by Equation (31). According to the analytical paradigm of spatial econometrics, the beta convergence trend cannot be determined by the estimated parameter of ${\widehat{\omega }}_1$. However, the beta convergence trend should be judged as the following three steps. First, we determine the data-generating process of the model of (30) and (31) and calculate the parameter marginal effect matrix, $S_1\left(W\right)=\frac{\partial Ln\left(GTF{P}_{i,t+1}/GTF{P}_{i,t}\right)}{\partial Ln\left(GTF{P}_{i,t}\right)}={\widehat{\Theta }}_2\left(I_{ST{W}_1}{\widehat{\omega }}_1+ST{W}_1{\widehat{\omega }}_2\right)$, wherein ^ indicates estimated values of corresponding parameters, ${\widehat{\Theta }}_2={\left(I_{ST{W}_1}-{\widehat{\omega }}_{0}ST{W}_1\right)}^{-1}$, and $I_{ST{W}_1}$ is an identity matrix with the same row and column as those of $ST{W}_1$. Second, we calculate the total effects, the direct effects, and the indirect effects based on the matrix of $S_1\left(W\right)$, wherein the total effects are defined as $\frac{1}{Nu{m}_{ST{W}_1}}{{\iota }^{\prime }}_{ST{W}_1}{S}_1\left(W\right){\iota }_{ST{W}_1}$, the direct effects are defined as $\frac{1}{Nu{m}_{ST{W}_1}}Trace\left[S_1\left(W\right)\right]$, and the indirect effects are defined as the remaining after subtracting the direct effects from the total effects. Moreover, $Trace\left(·\right)$ is the trace statistic, $Nu{m}_{ST{W}_1}$ is the dimensions of the spatiotemporal weight matrix, and ${\iota }_{ST{W}_1}$ indicates the vector with $Nu{m}_{ST{W}_1}$ rows and one column and its elements are always ones. Subsequently, the spatial convergence can be judged by the total effects. It shows an absolute beta convergence if the total effects are smaller than 0.

5.2.2. The Spatial Convergence Patterns of China’s Urban Industrial Green TFPs

According to Equation (29), sigma convergence is calculated for all 280 cities and for cities belonging to different spatial groups, including eastern, central, western, and northeast China, as depicted in Figure 8.

Overall, for all 280 cities, the coefficient of variation of China’s urban industrial green TFPs displayed dynamic adjustments before 2011, albeit showing an overall upward trend. However, after 2011, this coefficient of variation continuously declined. Overall, from 2003 to 2019, the coefficient of variation of China’s urban industrial green TFPs saw a slight decrease.

Looking at the sub-levels of different spatial groups, the coefficient of variation of industrial green TFPs in cities in eastern China exhibited noticeable fluctuations, with a slightly smaller coefficient of variation observed in 2019 compared to 2003. In contrast, the coefficient of variation of industrial green TFPs in cities in central China developed steadily and showed no significant changes during 2003–2019. Similarly, the coefficient of variation of industrial green TFPs in cities in western China mirrored the overall trend for all 280 cities, initially increasing before 2011 and then declining thereafter. The coefficient of variation of industrial green TFPs in cities in northeast China also followed a similar trend of initially increasing and then declining, albeit with the boundary year shifting to 2016.

Overall, based on the coefficients of variation, sigma convergence of China’s industrial green TFPs in the northeastern and central regions is not particularly evident. However, since 2011, sigma convergence of China’s industrial green TFPs has become apparent for all cities and for cities in western China.

Based on Equations (30) and (31), we can analyze the overall absolute beta convergence of China’s urban industrial green TFPs. As shown in Figure 9, the relationships of the variable’s pairs are figured, including $Ln\left(GTF{P}_{i,t+1}/GTF{P}_{i,t}\right)$ and $Ln\left(GTF{P}_{i,t}\right)$, $Ln\left(GTF{P}_{i,t+1}/GTF{P}_{i,t}\right)$, and both the spatial spillover terms of $ST{W}_1\times Ln\left(GTF{P}_{i,t+1}/GTF{P}_{i,t}\right)$ and $ST{W}_1\times Ln\left(GTF{P}_{i,t}\right)$. From Figure 9, there appear to be no discernible relationships between $Ln\left(GTF{P}_{i,t+1}/GTF{P}_{i,t}\right)$ and the three potential influencing factors. Consequently, this study proceeded with the trial estimation of the model outlined in Equations (30) and (31) and their degradation models, as detailed in

Table 6. Estimated results of the beta astringency models and their decomposed marginal effects.

		NSM	SXL	SAR	SEM	SDM	SDEM	SAC	GNSM
Const.		0.0473 (10.94) ***	0.0356 (4.78) ***	0.0032 (0.76)	0.0301 (2.91) ***	0.0127 (1.74) *	0.0122 (0.85)	0.0150 (2.13) **	0.0077 (0.56)
$Ln\left(GTF{P}_{i,t}\right)$		−0.0539 (−11.59) ***	−0.0495 (−9.57) ***	−0.0505 (−11.04) ***	−0.0562 (−11.14) ***	−0.0545 (−10.70) ***	−0.0571 (−10.01) ***	−0.0544 (−11.11) ***	−0.0555 (−10.95) ***
$ST{W}_1\times Ln\left(GTF{P}_{i,t}\right)$			−0.0204 (−1.93) *			0.0188 (1.81) *	−0.0414 (−2.07) **		−0.0146 (−0.69)
${\omega }_0$				0.4850 (121.70) ***		0.4990 (123.63) ***		0.3250 (6.87) ***	0.3000 (1.25)
${\omega }_5$					0.6180 (12.10) ***		0.6640 (5.24) ***	0.3510 (6.84) ***	0.5250 (2.79) ***
${\widehat{R}}^2$		0.0289	0.0295	0.0279	0.0619	0.0281	0.0642	0.0634	0.0657
${\sigma }^2$		0.0646	0.0646	0.0627	0.0624	0.0626	0.0622	0.0623	0.0621
$Log\left(L\right)$		−220.52	−218.65	1384.1	1397.52	1385.64	1401.32	1397.56	−144.6
the Beta astringency coefficients	Total effects	−0.0539	−0.0699	−0.0981	−0.0713	−0.0562	−0.0985	−0.0806	−0.1001
	Direct effects	−0.0539	−0.0495	−0.0507	−0.0546	−0.0562	−0.0571	−0.0545	−0.0556
	Indirect effects	0.0000	−0.0204	−0.0474	−0.0167	0.0000	−0.0414	−0.0261	−0.0445

Note: the above outputs were collected based on MATLAB R2023a. ( ) represents T-statistic, ***, **, and * mean having passed the hypothesis test with a significance level of 1%, 5%, and 10%.

. The results revealed abnormal estimated outcomes in models including SAR, SEM, SDM, SDEM, and SAC, as their log-likelihood values were positive, rendering these five models suboptimal. Additionally, the GNSM model could not be considered optimal due to the lack of significance in some estimated parameters. Furthermore, although the goodness of fit in the NSM and SXL models was relatively modest, the statistical characteristics of these two models were comparatively better. Upon decomposing the total effects from the marginal effects matrix, including the beta astringency coefficient in both the NSM and SXL models, it was observed that China’s urban industrial green TFPs exhibited a fundamental trend of absolute convergence in the long run.

The sub-dimensional absolute beta convergence of China’s urban industrial green TFPs, which belong to different spatial groups including eastern, central, western, and northeast China, has also been analyzed. The estimated parameters and their statistical characteristics of the corresponding models are listed in

Table 7. Estimated results of the beta astringency model in different spatial groups and their decomposed marginal effects.

			NSM	SXL	SAR	SEM	SDM	SDEM	SAC	GNSM
Cities in Eastern China	$Ln\left(GTF{P}_{i,t}\right)$		−0.0349 (−4.25) ***	−0.0314 (−3.61) ***	−0.0333 (−4.15) ***	−0.0403 (−4.69) ***	−0.0394 (−4.64) ***	−0.0397 (−4.51) ***	−0.0386 (−4.55) ***	−0.0405 (−4.68) ***
	$ST{W}_1\times Ln\left(GTF{P}_{i,t}\right)$			−0.0227 (−1.22)			0.0405 (2.24) **	0.0104 (0.38)		0.0378 (1.10)
	${\omega }_0$				0.5090 (73.48) ***		0.5709 (78.99) ***		0.3220 (3.80) ***	0.4759 (0.51)
	${\widehat{R}}^2$		0.0127	0.0131	0.0135	0.0639	0.0015	0.0631	0.0648	0.0687
	$Log\left(L\right)$		75.88	76.63	560.34	567.41	562.35	567.49	565.94	111.71
	the Beta astringency coefficients	Total effects	−0.0349	−0.0541	−0.0678	−0.0403	0.0026	−0.0293	−0.0569	−0.0052
		Direct effects	−0.0349	−0.0314	−0.0334	−0.0403	−0.0392	−0.0397	−0.0387	−0.0404
		Indirect effects	0.0000	−0.0227	−0.0344	0.0000	0.0418	0.0104	−0.0183	0.0352
Cities in Central China	$Ln\left(GTF{P}_{i,t}\right)$		−0.0424 (−5.23) ***	−0.0462 (−4.92) ***	−0.0406 (−5.04) ***	−0.0455 (−5.21) ***	−0.0507 (−5.43) ***	−0.0487 (−5.31) ***	−0.0428 (−5.11) ***	−0.0510 (−5.44) ***
	$ST{W}_1\times Ln\left(GTF{P}_{i,t}\right)$			0.0152 (0.80)			0.0423 (2.18) **	0.0230 (1.07)		0.0461 (1.58)
	${\omega }_0$				0.2010 (62.52) ***		0.2500 (5.38) ***		0.1340 (50.94) ***	0.2691 (0.37)
	${\widehat{R}}^2$		0.0202	0.0199	0.0188	0.0352	0.0196	0.0353	0.0350	0.0386
	$Log\left(L\right)$		74.65	74.98	523.69	527.42	525.90	527.99	525.98	86.62
	the Beta astringency coefficients	Total effects	−0.0424	−0.0310	−0.0508	−0.0455	−0.0112	−0.0257	−0.0494	−0.0067
		Direct effects	−0.0424	−0.0462	−0.0406	−0.0455	−0.0506	−0.0487	−0.0428	−0.0509
		Indirect effects	0.0000	0.0152	−0.0102	0.0000	0.0394	0.0230	−0.0066	0.0442
Cities in Western China	$Ln\left(GTF{P}_{i,t}\right)$		−0.0721 (−8.14) ***	−0.0701 (−7.22) ***	−0.0707 (−8.03) ***	−0.0748 (−7.95) ***	−0.0722 (−7.48) ***	−0.0740 (−7.74) ***	−0.0730 (−8.02) ***	−0.0729 (−7.24) ***
	$ST{W}_1\times Ln\left(GTF{P}_{i,t}\right)$			−0.0103 (−0.50)			0.0081 (0.39)	−0.0320 (−1.33)		−0.0117 (−0.14)
	${\omega }_0$				0.3270 (4.18) ***		0.3320 (4.13) ***		0.2220 (48.59) ***	0.2000 (0.20)
	${\widehat{R}}^2$		0.0469	0.0464	0.0450	0.0603	0.0446	0.0614	0.0600	0.0605
	$Log\left(L\right)$		−233.83	−233.71	232.10	233.65	232.16	234.56	233.47	−225.28
	the Beta astringency coefficients	Total effects	−0.0721	−0.0804	−0.1051	−0.0748	−0.0960	−0.1060	−0.0938	−0.1058
		Direct effects	−0.0721	−0.0701	−0.0708	−0.0748	−0.0723	−0.0740	−0.0731	−0.0730
		Indirect effects	0.0000	−0.0103	−0.0342	0.0000	−0.0237	−0.0320	−0.0208	−0.0328
Cities in Northeast China	$Ln\left(GTF{P}_{i,t}\right)$		−0.0934 (−5.92) ***	−0.0803 (−4.88) ***	−0.0880 (−5.61) ***	−0.0850 (−5.06) ***	−0.0795 (−4.87) ***	−0.0860 (−5.24) ***	−0.0869 (−5.47) ***	−0.0842 (−4.72) ***
	$ST{W}_1\times Ln\left(GTF{P}_{i,t}\right)$			−0.0841 (−2.64) ***			−0.0624 (−1.94) *	−0.1295 (−3.12) ***		−0.1176 (−1.73) *
	${\omega }_0$				0.3290 (3.34) ***		0.2540 (2.95) ***		0.2869 (3.13) ***	0.0716 (0.17)
	${\widehat{R}}^2$		0.0590	0.0693	0.0658	0.0704	0.0662	0.0878	0.0770	0.0871
	$Log\left(L\right)$		−60.74	−57.25	131.55	130.32	133.30	135.29	131.73	−53.19
	the Beta astringency coefficients	Total effects	−0.0934	−0.1644	−0.1311	−0.0850	−0.1902	−0.2155	−0.1219	−0.2174
		Direct effects	−0.0934	−0.0803	−0.0882	−0.0850	−0.0800	−0.0860	−0.0871	−0.0844
		Indirect effects	0.0000	−0.0841	−0.0429	0.0000	−0.1102	−0.1295	−0.0348	−0.1330

Note: the above outputs were generated using MATLAB R2023a. () represents T-statistic, ***, **, and * mean having passed the hypothesis test with a significance level of 1%, 5%, and 10%.

. Notably, in the process of total effects decomposition, the marginal effects matrices are determined solely by the spatiotemporal weight matrices, the spatial correlation coefficient of ${\widehat{\omega }}_0$, and the exogenous parameters of ${\widehat{\omega }}_1$ and ${\widehat{\omega }}_2$; thus, the estimated results of other variables are omitted in Table 7.

In Table 7, although the total effects of the beta astringency coefficient are negative for the cities in eastern China, there is still no conclusive evidence to prove the absolute convergence trend of China’s urban industrial green TFPs in the eastern region. This is due to two main reasons: first, the values of the log-likelihood of the eight potential models are all positive; second, the goodness of fit of the eight potential models is relatively low. Similar issues arise in the cities in central China, where there is also no evidence to confirm the absolute convergence trend of urban industrial green TFPs. For the cities in western China, the NSM is identified as the optimal model to assess absolute beta convergence because all parameters in the NSM are significant. In the NSM model of cities in western China, the value of log-likelihood is negative, and the total effects are −0.0349, indicating a certain beta convergence trend in industrial green TFPs. However, the goodness of fit for the NSM in western cities remains low, suggesting a weak beta convergence trend.

Regarding the cities in northeast China, both the NSM and the SXL can be employed to assess the beta convergence trend. There is indeed a beta convergence trend for industrial green TFPs in northeast cities, as indicated by the total effects of −0.0934 and −0.1644 for the NSM and SXL models, respectively. Moreover, the statistical characteristics of both the NSM and the SXL support the beta convergence trend. However, the beta convergence trend in northeast cities is also weak, given the small adjusted goodness of fit in the NSM and SXL models.

6. Conclusions and Comments

This paper introduces a novel spatiotemporal econometric approach called the Spatiotemporal Econometric Solow Residual Method (STE-SRM), which incorporates undesired outputs and spatial spillover terms to compute green Total Factor Productivities (TFPs). It applies this method to calculate the industrial green TFPs of 280 cities in China from 2003 to 2019. Additionally, it re-calculates China’s urban industrial green TFPs using both the DEA-SBM and the Bayesian SFA and evaluates the accuracy of the STE-SRM by comparing the results. Furthermore, the paper examines the spatial heterogeneity and spatial convergence of China’s urban industrial green TFPs using Dagum’s Gini coefficient, sigma astringency, and beta astringency analysis. The following conclusions are drawn:

(1) The STE-SRM enhances the traditional Solow residual method by integrating undesired outputs and spatial spillover terms. During construction, the sum of desired and undesired output shares is fixed at 1, and various spatial spillover terms are considered. While the calculation of industrial green TFPs can still follow the traditional Solow residual method under STE-SRM, determining input per capita shares and undesired outputs requires estimation and selection of empirical production function models.

(2) Utilizing the STE-SRM, the paper identifies the spatial Durbin error model with mixed effects as the optimal model for calculating China’s urban industrial green TFPs. Results show that industrial green TFPs in cities like Huangshan, Fangchenggang, and Sanya are relatively high, while those in Jincheng, Datong, and Taiyuan are comparatively low. Various factors such as industrial output per capita, industrial inputs per capita (including capital and energy), and ratios of undesired outputs (e.g., wastewater, sulfur dioxide, and smoke emissions) to desired outputs contribute to these discrepancies.

(3) Comparative analysis reveals that China’s urban industrial green TFPs computed by the STE-SRM exhibit a broader range of values and higher concentrations compared to those from the Bayesian SFA and DEA-SBM. Thus, the STE-SRM is recommended for researchers seeking wider range but concentrated industrial green TFPs, while the Bayesian SFA is suitable for those preferring narrower range but dispersed results.

(4) Spatial heterogeneity remains significant in China’s urban industrial green TFPs under the STE-SRM, with Dagum’s Gini coefficient ranging from 0.4 to 0.5. Despite similarities between intra-group and inter-group differences, differences in intensity of trans-variation contribute most to the coefficient. Since 2011, China’s urban industrial green TFPs have shown a clear sigma convergence trend, particularly evident in western cities. However, convergence trends in cities from the eastern, central, and northeastern regions are less apparent. Additionally, a slight beta convergence trend is observed in urban industrial green TFPs of western and northeastern cities, while such trends are less evident in cities from eastern and central China.

In this study, we introduce innovative methods for computing industrial green TFPs while maintaining the logic of the traditional Solow residual method. However, there are limitations that warrant further investigation. Firstly, while the STE-SRM improves upon traditional methods by incorporating spatial spillover terms and undesired outputs, future research should explore local spatial econometric approaches like geographical weighted regression (GWR) and geographical and temporal weighted regression (GTWR). Secondly, due to data limitations and inconsistencies, conditional beta convergence properties and the main factors affecting urban industrial green TFPs were not analyzed. These aspects require ongoing attention in subsequent research endeavors. Furthermore, this article delves into the comparison and interpretation of the precision of the STE-SRM, DEA-SBM, and Bayesian SFA methods. It primarily relies on data from the industrial sectors of 280 cities in China spanning from 2003 to 2019. It is worth noting that any fluctuations in the data could potentially impact the precision of these three methods, thus highlighting a deficiency in the theoretical underpinning of method precision interpretation within this article. Regrettably, the disparate theoretical foundations upon which these three methods are built contribute to a lack of a cohesive theoretical framework for analyzing and comparing measurement precision. This deficiency underscores a notable limitation of the article and serves as a pointer towards avenues for future research.