A Study of Total-Factor Energy Efficiency for Regional Sustainable Development in China: An Application of Bootstrapped DEA and Clustering Approach

Abstract

Total-factor energy efficiency (TFEE) is widely used to measure energy efficiency under the data envelopment analysis (DEA) framework, but the efficiencies obtained from DEA are structurally biased upward, and thus TFEE tends to overestimate energy efficiency. This research thus applies the bootstrapped DEA approach to correct the bias of TFEE. Using a dataset consisting of 30 provinces of China in the period 2016–2019, the bootstrapped-based test supports technology with variable returns to scale. The biased-corrected TFEE also indicates that energy consumption on average can be scaled down by 42.36%, rather than the biased value of 19.4%. The bootstrapped clustering partitions provinces into three groups: Cluster 1, with Guizhou as the representative medoid, includes half of the superior coastal provinces in terms of actual energy consumption and TFEE and half of the competitive inland provinces, whereas Cluster 3 outperforms Cluster 2 in terms of TFEE, but the actual energy consumption is higher, with Shandong and Hebei as the representative medoids, respectively. Lastly, empirical results imply that the northeast and central regions need more government attention and resources to practice sustainable development and improve TFEE.

1. Introduction

Since its reform and opening up in 1979, the People’s Republic of China (hereafter denoted as China) has achieved incredible development and is currently the world’s largest emerging economy and even the second-largest economy overall. Nevertheless, this rapid economic growth has led to excessive energy consumption and ecological problems [1,2]. From 1980 to 2014, its energy use per capita increased from 609.46 to 2224.36 kilotons of oil equivalent [3]. At the same time, China has also suffered severe air pollution [4], as its total greenhouse gas emissions rose from 2.76 to 12.34 million kilotons of CO₂ equivalent during the period from 1980 to 2018 [3].

The main idea of sustainable development is not only to meet the needs of the present without sacrificing the ability of future generations but also to achieve balanced development between human beings and nature [5]. In order to follow the trend of global ecological environment (eco-environment) protection and solve increasingly serious domestic environmental problems, China has actively promoted sustainable development in recent years. However, achieving a balance between regional economic growth and eco-environment protection is clearly a serious challenge for the country. Hence, the measurement of energy efficiency is a vital theme in the regional sustainable development of China, as the largest emerging economy in the world.

Eco-efficiency (ecological efficiency), which focuses on the efficient use of resources for economic and environmental advancement, is widely used to examine sustainability performance at industry, regional, and national levels [6,7,8]. Although it has attracted considerable attention from scholars [9,10], many studies ignore the possibility of multiple input substitution by using partial factor frameworks that lead to biased estimates [11]. For instance, the traditional energy efficiency index, which is one of the important eco-efficiency indicators, only regards energy as a single input and ignores the substitution or complementation of energy with other inputs such as labor and capital [12]. Thus, using a partial factor energy efficiency indicator may obtain a plausible result. Data envelopment analysis (DEA), a total-factor analysis approach, is essentially a linear programming model for performance evaluation by calculating the best multiplier for inputs and outputs. It has been widely applied in eco-efficiency assessments of regional sustainable development [11,13,14]. Hu and Wang [15] are the first to propose total-factor energy efficiency (TFEE) to calculate energy efficiency under the DEA framework.

DEA efficiencies, by construction, are biased upward [16], and thus TFEE inclines to misjudge energy efficiency. To our knowledge, no studies have used an appropriate method to correct for bias in TFEE. Furthermore, DEA models, according to the nature of production technology, can be generally classified into two categories: the CCR model (constant returns to scale, CRS) and the BCC model (variable returns to scale, VRS). When technology exhibits CRS, the eco-efficiency obtained from either the CCR model or the BCC mode is consistent, while the former is more efficient than the latter. Nevertheless, if the production technology is VRS, then the latter still remains consistent, but the former becomes inconsistent [17].

Conventional DEA approaches rely on linear programming techniques for solutions and thus cannot directly give clear statistical inferences. To properly evaluate the eco-efficiency and/or TFEE of regional sustainable development, we extend the bootstrapping procedure, initially proposed by Simar and Wilson [17], to statistically test the property of returns to scale for China’s regional sustainable development by considering quasi-fixed inputs and undesirable environmental emissions. After recognizing the appropriate returns to scale, we generate the appropriate bootstrapped samples to correct the bias of TFEE and then apply the bootstrapped clustering approach to classify provinces in the same group that are more similar to each other than to those in other groups.

The rest of the paper is organized as follows. Section 2 surveys the relevant literature. Section 3 explains the methodology. Section 4 provides a description of the data and empirical analyses. Section 5 presents discussion and policy implications. The final section is the concluding remarks.

2. Literature Review

Sustainable development must meet the needs of the present without sacrificing the ability of future generations and strive to achieve balanced development between human beings and nature. The concept of eco-efficiency provides a quantitative direction for strengthening intensive development while achieving sustainable development. The ratio between economic value and environmental load is the simplest way to evaluate eco-efficiency [18], but this single-factor framework is generally unsatisfactory in assessing the eco-efficiency of regional sustainable development [11]. For example, conventional energy efficiency, one of the important eco-efficiency indicators, only takes account of a single input such as energy and neglects the substitution or complement among energy, labor, and/or capital [11,12,19]. In other words, the partial factor energy efficiency indicator may deliver a plausible result. DEA, initially proposed by Charnes et al. [20] and extended by Banker et al. [21], is a linear programming model to evaluate efficiencies by calculating the optimal input and output multipliers. Because it can deal with multi-input and/or multi-output without specifying any functional form, it has been widely applied in measuring eco-efficiency of regional sustainable development [22,23].

DEA models essentially allow all inputs to be free and instantly adjustable to their optimum levels, thus ignoring the presence of non-radial slacks and resulting in weaker efficiency measures [24]. In view of this, Ali and Seiford [25] proposed a second-stage linear programming problem to maximize the sum of non-radial slacks. Hence, the target energy consumption is equal to the actual energy consumption minus the sum of the radial and non-radial energy consumptions. Hu and Wang [15] set up total factor energy efficiency (TFEE), which is the ratio of target energy consumption to actual input consumption, in order to measure energy efficiency under the DEA framework. Since then, many studies have chosen DEA models to measure TFEE [12,26,27,28].

Ouellette and Vierstraete [29] argued in the real world that not all inputs could be adjusted to their target levels due to regulation, adjustment costs, indivisibilities, etc. From the viewpoints of policymakers or the government, for instance, it is unacceptable to reduce the employment level and shrink knowledge stocks while reducing excess energy consumption. The quasi-fixed input DEA model handles both irreducible quasi-fixed inputs and reducible variable inputs simultaneously [29,30,31]. Shi et al. [31] treated non-energy inputs as quasi-fixed to measure regional industrial energy efficiency in China. Zhou and Ang [32] also conducted an energy efficiency index under a quasi-fixed input assumption.

In the actual production process, desirable outputs are generally accompanied by a range of undesirable or unwelcome by-products, yet traditional DEA models only deal with desirable outputs. Some studies used inverse transformations to convert the undesired output to the desired value [33] or applied the property of output translation invariance [34] to transform the negative values of undesirable outputs into positive values [35]. Other scholars assumed undesirable outputs to be weakly disposable [36,37,38].

DEA efficiencies are structurally biased upward [16] and thus tend to overestimate energy efficiencies. Some studies measure TFEE by considering quasi-fixed inputs and/or undesirable environmental emissions [39], but to our knowledge, there are no studies on regional sustainability using an appropriate method to correct the bias of TFEE. Therefore, extending the research on TFEE in China’s regional sustainable development via a bootstrapped DEA approach to measure biased-corrected TFEE is of great significance.

3. Methodology

Eco-efficiency focuses on the practices of efficient use of resources for economic and environmental progress, while industrial operations produce undesirable products such as sulfur dioxide (SO₂) emissions, which erode the environment. In addition, efficient use of labor input means increased unemployment, which is unwanted for policymakers and governments alike. Hence, an empirical model that can appropriately evaluate the eco-efficiency and/or TFEE of regional sustainable development should accommodate both undesirable outputs (such as SO₂ emissions) and irreducible quasi-fixed inputs (such as labor input).

Suppose that there are N DMUs. Each DMU employs k variable inputs $\underset{˜}{\mathit{x}}$ = $\left(x_1,\dots ,x_k\right)\in {\Re }_{+}^k$ and r quasi-fixed inputs $\underset{˜}{\mathit{v}}$ = $\left(v_1,\dots ,v_r\right)\in {\Re }_{+}^r$ to produce h undesirable outputs $\underset{˜}{\mathit{b}}$ = $\left(b_1,\dots ,b_h\right)\in {\Re }_{+}^h$ and m desirable outputs $\underset{˜}{\mathit{u}}$ = $\left(u_1,\dots ,u_m\right)\in {\Re }_{+}^m$. Charnes et al. [20] proposed the CCR model, using the smallest free disposal convex set covering all the data, in order to construct the production possibility set. Since the quasi-fixed inputs cannot be altered, the corresponding input-oriented eco-efficiency for the nth DMU is given by:

$${\widehat{\theta }}_n=\underset{\theta , \underset{˜}{\mathit{\lambda }}}{\inf }\{\theta |{\underset{˜}{\mathit{v}}}_n\ge V\underset{˜}{\mathit{\lambda }}, \theta {\underset{˜}{\mathit{x}}}_n\ge X\underset{˜}{\mathit{\lambda }}, -{\underset{˜}{\mathit{u}}}_n\ge U\underset{˜}{\mathit{\lambda }}, {\underset{˜}{\mathit{b}}}_n=B\underset{˜}{\mathit{\lambda }}, \underset{˜}{\mathit{\lambda }}\ge \underset{˜}{\mathbf{0}}\}$$

(1)

where $\theta$ is a scalar, $V=\left({\underset{˜}{\mathit{v}}}_1,\dots ,{\underset{˜}{\mathit{v}}}_N\right)$, $X=\left({\underset{˜}{\mathit{x}}}_1,\dots ,{\underset{˜}{\mathit{x}}}_N\right)$, $U=\left({\underset{˜}{\mathit{u}}}_1,\dots ,{\underset{˜}{\mathit{u}}}_N\right)$, $B=\left({\underset{˜}{\mathit{b}}}_1,\dots ,{\underset{˜}{\mathit{b}}}_N\right)$, $\underset{˜}{\mathit{\lambda }}$ is an (N × 1) vector of intensity variables, and $\underset{˜}{\mathbf{0}}$ is an (N × 1) vector of zeros. Note that the inequality constraints for $\underset{˜}{\mathit{x}}$, $\underset{˜}{\mathit{v}}$, and $\underset{˜}{\mathit{u}}$ reflect the characteristics of strong disposability, whereas the equality constraint for $\underset{˜}{\mathit{b}}$ reveals the property of weak disposability.

The CCR model assumes that the production exhibits constant returns to scale (CRS), which is only appropriate when all DMUs are operating at an optimal scale. Banker et al. [21] extended the CCR model to account for variable returns to scale (VRS), called the BCC model. Mathematically, the BCC model is modified easily from the CCR model by adding the convexity constraint ${\underset{˜}{\mathbf{1}}}^{\prime }\underset{˜}{\mathit{\lambda }}$ = 1 in Equation (1), where $\underset{˜}{\mathbf{1}}$ is an (N × 1) vector of ones.

If the production technology displays CRS globally, then both eco-efficiencies obtained from the CCR model, ${\widehat{\theta }}_n^c$, and the BCC model, ${\widehat{\theta }}_n^v$, are consistent, but the former is more efficient than the latter. Nevertheless, if the production technology exhibits VRS in some locations, then ${\widehat{\theta }}_n^v$ is still consistent, but ${\widehat{\theta }}_n^c$ is not [17]. Conventional DEA models use linear programming techniques to evaluate eco-efficiencies and thus cannot directly perform any statistical inference to investigate the property of returns to scale. Hence, this study extends the bootstrapping procedure of Simar and Wilson [17,40] by considering quasi-fixed inputs and undesirable outputs to generate bootstrap samples, thereby testing the returns to scale. We summarize the bootstrapping procedure as follows.

(1)

Use all DMUs to evaluate ${\widehat{z}}_n$ (=$1/{\widehat{\theta }}_n$), n = 1,…, N, by an appropriate DEA model.

(2)

Randomly draw a sample of size N with replacement from the following set:

$$\left\{\left(2-{\widehat{z}}_1\right),\dots ,\left(2-{\widehat{z}}_N\right),{\widehat{z}}_1,\dots , {\widehat{z}}_N\right\}$$

To obtain $\{z_1^{*},\dots , z_N^{*}\}$.

(3)

Compute ${\tilde{z}}_n^{*}={\overline{z}}^{*}+\left({\delta }_n^{*}-{\overline{z}}^{*}\right)/{[1+{\left(\tau /{\widehat{\sigma }}^{*}\right)}^2]}^{0.5}\begin{array}{cc},& n=1,\dots ,N\end{array}$, where ${\overline{z}}^{*}$ and ${\widehat{\sigma }}^{*}$ are respectively the mean and standard deviation of $\{z_1^{*},\dots , z_N^{*}\}$, ${\delta }_n^{*}=z_n^{*}+\tau {\omega }_n$, ${\omega }_n$ draws randomly from a standard normal distribution, $\tau =1.06min\{\widehat{\sigma }, \widehat{R}/1.34\}{N}^{-1/5}$, and $\widehat{\sigma }$ and $\widehat{R}$ are respectively the standard deviation and the interquartile range of $\left\{{\widehat{z}}_1,\dots , {\widehat{z}}_N\right\}$.

(4)

Set ${\underset{˜}{\mathit{x}}}_n^{*}={\tilde{z}}_n^{**}{\underset{˜}{\mathit{x}}}_n/{\widehat{z}}_n$, where ${\tilde{z}}_n^{**}=2-{\tilde{z}}_n^{*}$ if ${\tilde{z}}_n^{*}<1$ and ${\tilde{z}}_n^{**}={\tilde{z}}_n^{*}$ otherwise.

(5)

Compute the bootstrap estimate ${\widehat{\theta }}_n^{*}$ (n = 1,…, N) by the appropriate DEA model with technology (V, X*, U, B), where X*$=\left({\underset{˜}{\mathit{x}}}_1^{*},\dots ,{\underset{˜}{\mathit{x}}}_N^{*}\right)$.

(6)

Re-do steps (2)–(5) T times to achieve bootstrap estimates ${\{{\widehat{\theta }}_{nt}^{*}\}}_{n=1}^N$ for t = 1,…, T.

Since the assumption of CRS is more restrictive than that of VRS, we follow Simar and Wilson’s [40] suggestion by setting CRS as the null hypothesis H₀ and the test statistic as:

$$\mathsf{\Pi }=N^{-1}{\displaystyle \sum }_{n=1}^N{\widehat{\theta }}_n^c/{\widehat{\theta }}_n^v$$

(2)

To perform the bootstrapping-based test of returns to scale, we first calculate eco-efficiencies in step (1) by the CCR model and then compute the bootstrap estimate ${\mathsf{\Pi }}_t^{*}$ for each bootstrap sample t to obtain the bootstrap sample ${\left\{{\mathsf{\Pi }}_t^{*}\right\}}_{t=1}^T$. Because ${\widehat{\theta }}_n^c/{\widehat{\theta }}_n^v\le 1$ by construction and ${\widehat{\theta }}_n^c\approx {\widehat{\theta }}_n^v$ under H₀, we will reject H₀ if the test statistic $\mathsf{\Pi }$ is small enough. The p-value can be approximated by the proportion of bootstrap samples with values ${\mathsf{\Pi }}_t^{*}$ less than ${\pi }_0$ (the observed value of $\mathsf{\Pi }$ under H₀ and the observed dataset):

$$p\text{-}\mathrm{value}\approx {{\displaystyle \sum }}_{t=1}^T\mathrm{I}\left({\mathsf{\Pi }}_t^{*}\le {\pi }_0\right)/T,$$

(3)

where I(η) = 1 if η is true and 0 otherwise.

The total-factor energy efficiency (TFEE), initially proposed by Hu and Wang [15], is widely used to calculate the relative efficiency index of energy usage under the DEA framework. By definition, TFEE is a ratio of target energy input to actual energy input:

$$\mathrm{TFEE}=\frac{\mathrm{Target} \mathrm{energy} \mathrm{input}}{\mathrm{Actual} \mathrm{energy} \mathrm{input}}$$

(4)

Because target energy input is always less than or equal to actual energy input, the range of TFEE is between 0 and 1.

The DEA model assumes that all inputs are radially (proportionally) adjustable and tries to scale down all inputs proportionally as much as possible. However, redial input slacks, the difference between actual inputs and the redial target inputs (actual input minus $\widehat{\theta }$* actual input), may not always identify all (efficient) input slacks. Ali and Seiford [25] proposed a second-stage linear programming problem to find all possible non-radial input and output slacks. Therefore, the target input is equal to the actual input minus the sum of the radial and non-radial inputs. In other words, the total-factor energy efficiency is:

$$\mathrm{TFEE}=1-\frac{ \mathrm{Radial} \mathrm{energy} \mathrm{slack}+\mathrm{Non}-\mathrm{radial} \mathrm{energy} \mathrm{slack}}{\mathrm{Actual} \mathrm{energy} \mathrm{input}}.$$

(5)

4. Empirical Analysis

4.1. Data Sources and Input-Output Variables

The dataset, obtained from China Statistical Yearbook, the China Energy Statistical Yearbook, and China Statistical Yearbook on Environment, consists of 30 provinces (12 from the coastal region and 18 from the inland region) and 120 observations in China for the period 2016–2019. We exclude Tibet because of missing data. Since we have 4 years of data, all nominal variables are deflated by the GDP deflator with 2015 as the base year in order to remove price variation.

Choosing appropriate input and output variables is one of the key tasks in applying DEA models. To measure energy consumption, we first multiply the amount of each energy by the corresponding standard coal reference coefficient and then add it up to achieve the energy consumption per ton of standard coal equivalent (Energy). The data come from the China Energy Statistical Yearbook.

Classical production theory states that both labor and capital are key inputs for producing goods and services [41,42]. Due to data unavailability on capital stock (Capital), we apply the perpetual inventory method to measure capital stock $K_t$ in period t by:

$$K_t=\frac{\mathrm{Gross} \mathrm{capital} \mathrm{formation} \mathrm{in} \mathrm{period} t}{{\beta }_t+{\delta }_K},$$

where ${\delta }_K$ is the depreciation rate of capital stock, assumed to be 0.05 [43], and ${\beta }_{t }$ is the growth rate of output in period t. Since policymakers are generally reluctant to decrease employment levels and shrink capital stock while reducing excess energy consumption [31,32,39], we treat the number of employees (labor) and capital stock (capital) as irreducible quasi-fixed inputs.

Endogenous growth theory [44,45] encourages the government and the private sector to nurture innovation initiatives and encourages incentives for individuals and enterprises to be more creative, such as research and development (R&D) funding. For sustainable development, R&D stock is obviously not only a key input and nor should it be considered to reduce its quantity. Hence, this study regards R&D stock (RD) as a quasi-fixed input. Because data on R&D stock are unavailable, we also use the perpetual inventory method to measure R&D stock $R{D}_t$ in period t by:

$$R{D}_t=\frac{\mathrm{R}\&\mathrm{D} \mathrm{expenditure} \mathrm{in} \mathrm{period} t}{{\beta }_t+{\delta }_{RD}},$$

where ${\delta }_{RD}$ is the depreciation rate of R&D stock, assumed to be 0.15 [46,47,48].

The desirable outputs consist of GDP (GDP) and patents filed (Patent) of various provinces. Sulfur dioxide (SO₂) emissions are a serious problem in air pollution, especially in industrial clusters [49], and excessive emissions severely restrict sustainable development in China [50]. Hence, we treat SO₂ emissions (SO₂) as an undesirable output [50,51]. In sum, this study includes one input (Energy), three quasi-fixed inputs (Labor, Capital, and RD), one undesirable output (SO₂), and two desirable outputs (GDP and Patent).

Table 1. Descriptive statistics of inputs and outputs.

Variables	Mean	Std. Dev.	Min	Max
Variable input
Energy (10 million tons)	7.867	4.482	0.703	20.331
Quasi-fixed inputs
Labor (million)	27.731	18.091	3.243	71.503
Capital (RMB trillion)	19.552	13.533	2.636	60.053
RD (RMB trillion)	0.203	0.258	0.003	1.210
Undesirable output
SO2 (10,000 tons)	20.315	14.316	0.192	72.976
Desirable outputs
GDP (RMB trillion)	3.144	2.468	0.261	11.931
Patent (1000)	125.691	160.752	3.181	807.700

Note: All nominal variables are deflated by the GDP deflator with 2015 as the base year.

reports the summary statistics of the inputs and outputs used in the analysis.

The output and input variables in the DEA model should satisfy the property of isotonicity; i.e., increased inputs cannot reduce outputs.

Table 2. Correlation coefficients between outputs and inputs.

	Energy	Capital	Labor	RD
GDP	0.6813 (<0.001)	0.8362 (<0.001)	0.8208 (<0.001)	0.9648 (<0.001)
Patent	0.4578 (<0.001)	0.6087 (<0.001)	0.6444 (<0.001)	0.9359 (<0.001)
SO2	0.6224 (<0.001)	0.2857 (0.002)	0.3608 (<0.001)	0.1563 (0.088)

Note: The values in parentheses are p-value.

presents the Pearson correlation coefficients between input and output variables. All values are indeed positive. Furthermore, they are significant at the 1% level, except for the value between SO₂ and RD with a p-value = 0.088. Hence, our selected input and output variables explicitly meet the property of isotonicity.

4.2. Bootstrapped-Based Test of Returns to Scale

A critical issue in the DEA model is whether the production technology has constant or variable returns to scale. If the production possibility set exhibits CRS globally, then both eco-efficiencies obtained from the CCR model, ${\widehat{\theta }}^c$, and the BCC model, ${\widehat{\theta }}^v$, are consistent estimators, but the former is more efficient than the latter. This implies that the CCR model is more suitable for evaluating eco-efficiency than the BCC model. On the other hand, if the underlying technology displays VRS at some locations, then ${\widehat{\theta }}^c$ turns out to be inconsistent, while ${\widehat{\theta }}^v$ still remains consistent [17]. This situation suggests that we should employ the BCC model to measure eco-efficiency. Since traditional DEA models use linear programming techniques to calculate efficiencies, they cannot generally provide a formal statistical test of returns to scale. We extend Simar and Wilson’s [17] bootstrap procedure by considering quasi-fixed inputs and undesirable environmental emissions to perform a statistical test of returns to scale.

This study uses R software for empirical analysis.

Table 3. Bootstrapped-based test of returns to scale.

H0: CRS. Ha: VRS	${\pi }_o$	Critical Values		p-Value
		α = 0.01	α = 0.05
	0.8701	0.8854	0.8947	0.0005

Note: We use 2000 replications for the bootstrapped procedure.

presents that based on 2000 replications, the observed test statistic (${\pi }_o$= 0.8701) is less than the critical value of 0.8854 for a 1% level test (p-value = 0.0005). Hence, we reject the CRS hypothesis and use the BCC model to measure eco-efficiency and thereby investigate the relevant properties of China’s regional sustainable development.

4.3. Analysis of Total-Factor Energy Efficiency

Total-factor energy efficiency (TFEE), proposed by Hu and Wang [15], is widely used to measure energy efficiency under a total-factor framework. By definition, TFEE is a ratio of target energy input to actual energy input. Since only the energy input can be adjusted and the other inputs are quasi-fixed, there is no non-radio energy slack, which means that TFEE and eco-efficiency are the same [39]—that is:

$$\mathrm{TFEE}=1-\frac{(1-\widehat{\theta })\ast \mathrm{Actual} \mathrm{energy} \mathrm{input}}{\mathrm{Actual} \mathrm{energy} \mathrm{input}}=\widehat{\theta }.$$

The eco-efficiency obtained from the CCR or BCC model is structurally biased upward [16], as is TFEE naturally. Hence, it is desirable to use the bootstrap process to correct this upward bias. Let ${\{{\widehat{\theta }}_{nt}^{*}\}}_{t=1}^T$ be the bootstrap sample for $n=1,\dots ,N$. The bias-corrected estimator ${\widehat{\theta }}_n^{bc}$ is:

$${\widehat{\theta }}_n^{bc}={\widehat{\theta }}_n-\widehat{\mathsf{\Delta }}({\widehat{\theta }}_n)$$

(6)

where $\widehat{\mathsf{\Delta }}({\widehat{\theta }}_n)={\overline{\theta }}_n^{*}-{\widehat{\theta }}_n$, and ${\overline{\theta }}_n^{*}$ is the average value of the bootstrap sample ${\{{\widehat{\theta }}_{nt}^{*}\}}_{t=1}^T$. However, Efron and Tibshirani [52] pointed out that correcting for the bias also introduces additional noise and proposed that the condition for bias correction is:

$$|\widehat{\mathsf{\Delta }}({\widehat{\theta }}_n)|>se({\widehat{\theta }}_n^{*})/4$$

(7)

where $se({\widehat{\theta }}_n^{*})={[{\displaystyle \sum }_{t=1}^T{({\widehat{\theta }}_{nt}^{*}-{\overline{\theta }}_n^{*})}^2/T]}^{1/2}$.

Following the suggestion of Simar and Wilson [40], this study uses the bootstrap procedure for 2000 replications to obtain the bootstrap sample and to correct the bias. Since bootstrapping results do fulfill the condition of Equation (7) for each DMU (the mean values of $|\widehat{\mathsf{\Delta }}\left({\widehat{\theta }}_n\right)|$ and $se({\widehat{\theta }}_n^{*})/4$ are 0.2296 and 0.1902, respectively), we use Equation (6) to correct the bias of TFEE.

Let ${\mathrm{TFEE}}_{qf}$ and ${\mathrm{TFEE}}_{qf}^{bc}$ be the total-factor energy efficiencies equal to the BCC model eco-efficiency (=$\widehat{\theta }$) and the bias-corrected eco-efficiency (=${\widehat{\theta }}^{bc}$), respectively.

Table 4. Descriptive statistics of total-factor energy efficiency.

	Mean	Median	Std. Dev.	Min	Max
${\mathrm{TFEE}}_{qf} (=\widehat{\theta })$	0.8060	0.8648	0.2073	0.3339	1
${\mathrm{TFEE}}_{qf}^{bc}$$(={\widehat{\theta }}^{bc})$	0.5764	0.5812	0.1031	0.2998	0.8299

Notes: (i) We use 2000 replications for the bootstrapped procedure. (ii) Bootstrapping results do fulfill the condition of Equation (7) for each DMU, and thus we correct the bias of eco-efficiencies.

presents the descriptive statistics of ${\mathrm{TFEE}}_{qf}$ and ${\mathrm{TFEE}}_{qf}^{bc}$. The average value of ${\mathrm{TFEE}}_{qf}$ is 0.806, implying on average that energy consumption could drop 19.4%. However, since ${\mathrm{TFEE}}_{qf}$ is biased upward, the mean value of ${\mathrm{TFEE}}_{qf}^{bc}$ at 0.5764 indicates on average that energy consumption can be scaled down by 42.36% instead of 19.4%. Furthermore, 42 DMUs (35%) are BCC-efficient ($\widehat{\theta }$ = 1), whereas the largest value of bias-corrected eco-efficiency ${\widehat{\theta }}^{bc}$ is 0.8299, indicating that there is no efficient DMU operating on the production frontier after correcting for bias.

4.4. Cluster Analysis

The aim of cluster analysis is to group similar provinces in the same cluster based on a set of variables. Several studies employed it to supplement the DEA literature [37,53,54]. The k-medoids or partitioning around medoids (PAM) algorithm is a clustering approach to map a distance matrix into k clusters. It is based on the search for k representative provinces, medoids, among the provinces of a dataset [55]. Unlike the k-means algorithm, the k-medoids clustering minimizes a sum of pairwise dissimilarities instead of a sum of squared Euclidean distances. Hence, it is more robust to noise and outliers than the k-means algorithm. Moreover, the medoids are robust representations of the cluster centers, while the cluster center in k-means is not necessarily one of the provinces (it is the average between the provinces in the cluster).

The bootstrap samples of eco-efficiency offer sufficient information to construct a distance matrix based on the likelihood that the rank of the two provinces may be reversed [37,56]. Let ${\overline{\theta }}_n$ be the mean eco-efficiency of the nth province over the sample periods. If the nth province absolutely outperforms the ith province, then Pr(${\overline{\theta }}_n-{\overline{\theta }}_i>0$) = 1. Contrarily, if the nth province performs completely worse than the ith province, then Pr(${\overline{\theta }}_n-{\overline{\theta }}_i>0$) = 0. Both situations show maximum dissimilarity between the nth province and the ith province because both provinces never change in rank over all of the bootstrap samples.

In contrast to the above two extreme situations, Pr(${\overline{\theta }}_n-{\overline{\theta }}_i>0$) = 0.5 reveals the minimum dissimilarity between the nth province and the ith province because both provinces have an equal probability of outperforming each other. Therefore, we define the distance between the nth province and the ith province by:

$$D_{ni}=200\ast |\mathrm{Pr}({\overline{\theta }}_n-{\overline{\theta }}_i>0)-0.5|$$

(8)

Higher values of $D_{ni}$ indicate a bigger dissimilarity between the nth province and the ith province. Let $D_{ni}$ be the scalar at row n and column i of the (symmetric) dissimilarity matrix ${\left[D_{ni}\right]}_{N\times N}$. The k-medoids clustering minimizes the average dissimilarity of provinces to their closest representative province (medoid).

Figure 1–c are bivariate cluster plots for k = 2, k = 3, and k = 4, respectively. All 30 provinces are represented by points in the plot, using the first two principal components. Around each cluster is drawn an ellipse. Figure 1a shows for k = 2 that Cluster 1 contains 1.5 times as many provinces as Cluster 2 (18 vs. 12). Furthermore, Figure 1b,c reveal for k = 3 and k = 4 that Cluster 1 contains the same elements, while the elements of Cluster 4 of Figure 1c come evenly from Cluster 2 and Cluster 3 in Figure 1b.

Table 5. Mean values of partitioning around medoids (PAM) for k = 2, 3, 4.

	k = 2		k = 3			k = 4
Cluster	1	2	1	2	3	1	2	3	4
${\mathrm{TFEE}}_{qf}$	0.927	0.625	0.977	0.618	0.830	0.977	0.585	0.862	0.750
${\mathrm{TFEE}}_{qf}^{bc}$	0.614	0.520	0.607	0.521	0.611	0.607	0.489	0.637	0.599
Energy	7.932	7.738	6.550	7.409	10.261	6.550	7.611	9.501	9.517

presents mean values of partitioning around medoids for k = 2, 3, 4. For k = 2, both Clusters have similar actual energy consumption levels (7.932 vs. 7.738), but Cluster 1 outperforms Cluster 2 (0.614 vs. 0.520). Columns 4 to 6 show that Cluster 3 has the highest actual energy consumption (10.261), and Cluster 1 and Cluster 2 maintain similar levels (6.550 vs. 7.409), while both Cluster 1 and Cluster 3 have similar ${\mathrm{TFEE}}_{qf}^{bc}$ levels (0.607 vs. 0.611) and outperform Cluster 2 (0.521). These facts may suggest for k = 3 that the elements of Cluster 3 are mainly from Cluster 1 (for k = 2) with a higher actual energy consumption level and Cluster 2 (for k = 2) with both higher TFEE and actual energy consumption levels. For k = 4, Cluster 3 and Cluster 4 are very similar. Therefore, we choose k = 3 to perform a cluster analysis on the TFEE of China’s regional sustainable development.

Table 6. Results of partitioning around medoids (PAM) for k = 3.

	Cluster 1			Cluster 2			Cluster 3
	${\mathbf{TFEE}}_{\mathit{q}\mathit{f}}$	${\mathbf{TFEE}}_{\mathit{q}\mathit{f}}^{\mathit{b}\mathit{c}}$	Energy	${\mathbf{TFEE}}_{\mathit{q}\mathit{f}}$	${\mathbf{TFEE}}_{\mathit{q}\mathit{f}}^{\mathit{b}\mathit{c}}$	Energy	${\mathbf{TFEE}}_{\mathit{q}\mathit{f}}$	${\mathbf{TFEE}}_{\mathit{q}\mathit{f}}^{\mathit{b}\mathit{c}}$	Energy
Mean	0.977	0.607	6.550	0.618	0.521	7.409	0.830	0.611	10.261
Min	0.920	0.582	0.770	0.452	0.390	3.429	0.707	0.515	5.067
Max	1.000	0.631	15.046	0.767	0.672	18.958	0.943	0.731	17.843
Representative province		Guizhou			Hebei		Shandong
	0.983	0.625	5.475	0.560	0.401	18.958	0.821	0.565	17.843
Number of provinces Coast (inland)		11 6 (5)			11 3 (8)			8 3 (5)

shows the results of partitioning around medoids for k = 3. Cluster 2 has the lowest mean values for both BCC TFEE (${\mathrm{TFEE}}_{qf}$ = 0.618) and bias-corrected TFEE (${\mathrm{TFEE}}_{qf}^{bc}$ = 0.521), while Cluster 3 owns the highest actual energy consumption (10.261), and its mean ${\mathrm{TFEE}}_{qf}^{bc}$ is close to Cluster 1 (0.611 vs. 0.607). In addition, Cluster 1 not only contains all BCC eco-efficient provinces (Jiangsu, Guangdong, Sichuan, and Qinghai) that were consistently operating at the BCC frontier ($\widehat{\theta }$ = 1) during the sample period, but also all TFEEs are greater than the mean and median, regardless of ${\mathrm{TFEE}}_{qf}$ or ${\mathrm{TFEE}}_{qf}^{bc}$. In summary, we may characterize the samples of Cluster 1 as half being superior coastal provinces in terms of actual energy consumption and TFEE and half of those competitive inland provinces. The rest are divided into Cluster 2 and Cluster 3—the latter outperforms the former in terms of TFEE, but with a higher level of actual energy consumption.

One useful practice of the k-medoids clustering applications is to characterize the clusters by means of typical or representative objects. Because medoids are the representative province of a cluster and also observable objects, they can be used conveniently and/or economically to represent that cluster—that is, these medoids propose a research opportunity by using only a small set of k medoids instead of a large set one started off with. Table 6 shows that the representative medoids for Cluster 1, Cluster 2, and Cluster 3 are Guizhou, Hebei, and Shandong, respectively. This implies that we only need to investigate these three provinces to appropriately capture the key characteristics of TFEE of China’s regional sustainable development.

5. Discussion and Policy Implications

The single-factor framework is generally insufficient to assess energy efficiency because of neglecting substitution or complementarity among energy and other inputs such as labor and capital [11,19]. For instance, the mean ratios of actual energy consumption to GDP for Cluster 1, Cluster 2, and Cluster 3 are 3.28, 2.99, and 3.57, respectively. Cluster 3 has the highest mean ratio, but its TFEE is similar to Cluster 1 and better than Cluster 2. Hence, the partial factor energy efficiency indicator may bring a specious result. The government should avoid using the results of single-factor frameworks to formulate energy and/or sustainability policies.

DEA-based TFEE tends to overestimate energy efficiency because DEA efficiencies are structurally biased upward. Empirical results indicate, on average, that China’s energy consumption could be reduced by 42.36% (${\mathrm{TFEE}}_{qf}^{bc}$ = 0.5764), rather than the biased value 19.4% (${\mathrm{TFEE}}_{qf}$ = 0.806). This may suggest that previous DEA-based TFEE studies [12,26,27,28] are very likely to overvalue energy efficiency, thereby underestimating the energy consumption that can be decreased. Therefore, authorities could consider increasing energy-saving efforts if their policies follow DEA-based TFEEs.

Figure 2 exhibits provinces in Cluster 1, Cluster 2, and Cluster 3. It appears that the TFEE of the Yangtze River Delta and southwest belt regions (the main areas covered by Cluster 1) is better than that of other regions. This may suggest a U-shape relation between TFEE and economic development in areas of China, as the Pearl River Delta (in Guangdong) and the Yangtze River Delta became highly developed regions earlier, while western provinces such as Xinjiang, Qinghai, and Guizhou are underdeveloped areas. In addition, the government should aggressively improve the energy efficiency of Cluster 3, mainly in the central and northeastern regions, including Gansu, Chongging, Shaaxi, Hubei, Hunan, and Anhui (central region), Jilin, and Heilongjiang (northeast region), and Fujian, Hebei, and Tianjin (coastal region).

6. Conclusions

China has experienced a high rate of economic growth since its reform and opening-up in 1978, as seen by GDP rising quickly from $0.19 to $14.28 trillion from 1980 to 2019. Nevertheless, this booming economy has paid its dues in massive resource consumption and environmental pollution, which have put brakes on social and economic development. Considering undesired emission output and quasi-fixed inputs, we employ the bootstrapping procedure to measure TFEE of China’s regional sustainable development. After recognizing correct returns to scale by the bootstrapped-based test, this study uses a bootstrapping approach to correct the bias of TFEE and perform cluster analysis.

The dataset consists of 30 provinces (excluding Tibet because of missing data) over the period 2016–2019, thus including 120 observations. The test for returns to scale supports production technology to be VRS. In addition, biased-corrected TFEE indicates that energy consumption can be reduced by 42.36% instead of 19.4%. The k-medoids clustering, based on the bootstrap samples of TFEE, partitions China’s provinces into three clusters: Cluster 1 includes half of the superior coastal provinces in terms of actual energy consumption and TFEE and half of the competitive inland provinces, while Cluster 2 performs worse than Cluster 3, but with lower actual energy consumption. Furthermore, the representative medoids of Cluster 1, Cluster 2, and Cluster 3 are Guizhou, Hebei, and Shandong, respectively. Empirical results suggest that the government should actively promote sustainable development and TFEE in the northeast and central regions.