# Industrial Efficiency Algorithm Based on Spatio-Temporal-Data-Driven

Authors: Hongqu Lv; Wensi Cheng
Journal: Wireless Communications and Mobile Computing
 (2021)

---

## Abstract

Stochastic frontier model is an important and effective method to calculate industry efficiency. However, when dealing with temporal and spatial data from the industry, it is difficult to accurately calculate the industrial production efficiency due to the influence of spatial correlation and time lag effect. If the traditional spatial statistical method is used, the setting method of spatial weight matrix is often questioned. To solve this series of problems, one possible idea is to design a spatial data mining process based on stochastic frontier analysis. Firstly, the stochastic frontier model should be improved to analyze spatio-temporal data. In order to accurately measure the technical efficiency in the case of dual correlation between time and space, a more effective spatio-temporal stochastic frontier model method is proposed. Meanwhile, based on the idea of generalized moment estimation, an estimation method of spatiotemporal stochastic frontier model is designed, and the consistency of estimators is proved. In order to ensure that the most appropriate spatial weight matrix can be selected in the process of model construction, theK-fold crossvalidation method is adopted to evaluate the prediction effect under the data-driven idea. This set of spatio-temporal data mining methods will be used to measure the technical efficiency of high-tech industries in various provinces of China.


---

## Body

## 1. Introduction

Stochastic frontier analysis (SFA) is an important method to measure technical efficiency and calculate total factor productivity. The whole process is divided into two steps: the first step is the model estimation process, which can be regarded as a supervised learning process; the second step is to use the estimated model to calculate the technical efficiency, which can be regarded as an unsupervised learning process.From the perspective of machine learning, supervised processes have three main objectives: (a) feature selection and reduction of the dimension of feature variables; (b) selecting the optimal one from multiple classifiers or prediction models; (c) model evaluation, which estimates the prediction error of the selected classifier or prediction model on the new data.The paper found that the traditional stochastic frontier analysis method has the following defects: (a) it is not suitable for the special structure of spatial data or spatio-temporal data; (b) the modeling process lacks variety. The traditional analysis process is knowledge-driven and completely relies on a single theoretical model for estimation and testing. The above two characteristics lead to the large deviation of the traditional stochastic frontier model when analyzing the spatio-temporal data, and it is impossible to make an accurate measure of the production efficiency with spatial relationship, either. To solve the two problems above, this study considers two improvements to the industrial efficiency calculation process based on temporal and spatial data: (1) improve the existing stochastic frontier model and make it suitable for spatial data or spatio-temporal data; (2) turn the modeling process into a spatial data mining process. In view of the unique structure of spatio-temporal data, a more suitable crossvalidation method is proposed for the selection of prediction model.Stochastic frontier analysis (SFA) was successively proposed by Aigner et al. (1977) [1], Meeusen and Broeck (1977) [2], and Battese and Corra (1977) [3]. Over the past 40 years, its theoretical system and methods have been continuously expanded and innovated; it is widely used to measure the operating efficiency of different industries.The development of spatial statistics provides a theoretical basis for studying spatial interactions in stochastic frontier models. Druska and Horrace (2004) [4] first applied the method of spatial econometrics to the analytical framework of stochastic frontier model and started the research of spatial stochastic frontier model. Affuso (2010) [5] established the spatial stochastic frontier model and gave the maximum likelihood estimation in the empirical study. Tonini and Pede (2011) [6] applied maximum entropy method to parameter estimation of spatial stochastic frontier model. Vidolia et al. (2016) [7], Tsionas and Michaelides (2016) [8], Carvalho (2018) [9], and Adetutu et al. (2015) [10] consider SF models with local spatial dependence. Jin and Lee (2020) [11] proved the asymptotic properties of a maximum likelihood estimator of a spatial autoregressive stochastic frontier model. Kutlu et al. (2020) [12] proposed a spatial autoregressive stochastic frontier model, which allows for the endogeneity in both the frontier and environmental variables, and discussed a single-stage control function approach to estimate the parameters.Because spatial stochastic frontier analysis methods can fully consider the impact of spatial correlation, they can obtain more accurate results in efficiency analysis of data with spatial spillover effect and thus have been more widely used in recent years. Bergantino et al. (2020) [13] analyses the potential impact of airport competition on technical efficiency by applying the spatial stochastic frontier. Graaff (2020) [14] used spatial stochastic frontier model to estimate spatially correlated technical efficiencies within a European regional production function context. At present, some literatures have studied panel spatial stochastic frontier model, for example, Druska and Horrace (2004) [4], Tonini and Pede (2011) [6], and Lin Jia-Xian (2014) [15]. These literatures all focus on the static panel space stochastic frontier model, and the model utilizes two-dimensional information from panel data; formally, the spatial lag term of the explained variable and the spatial lag term of the error are used to capture the spatial correlation of the production unit. The time lag term is not included in the model, which means that the model still cannot fit well when there is significant inertia in the research problem. In input-output analysis, current behavior is largely dependent on past behavior, for example, the adjustment of capital stock is often influenced by previous capital. Therefore, a dynamic stochastic frontier model should be established, and the model should describe the double lag effect of space and time, so as to reflect the influence relationship between economic variables more objectively. The spatial weight matrix in spatial statistics is often considered to be “subjective.” Moreover, due to the various setting methods of spatial weight matrix, the selection of different spatial weight matrix may lead to the difference of model estimation results. In addition, the selection of spatial weight matrix has not formed a unified principle. Based on the above three points, the spatial weight matrix is often questioned. But in the era of “big data,” such skepticism may end [16].This paper proposes the spatiotemporal stochastic frontier model; considering that the model may be endogenous in time and space dimensions, a generalized method of moments (GMM) estimation process is designed to estimate the model. When Druska and Horrace (2004) [4] studied the static panel space stochastic frontier model, a generalized moment estimation process was proposed by referring to Kelejian and Prucha (1999) [17] for spatial error correlation. In this paper, Druska and Horrace (2004) [4] is used to deal with model’s error space autocorrelation, which is different from that of the stochastic frontier model. According to the method of Kapoor et al. (2007) [18], the compound error term was processed, and the moment condition was constructed to estimate the distribution parameters of the error term. In this paper, Jacobs et al. (2009) [19] was used as a reference to construct the moment condition, and Anselin (1988) [20] was used as a reference for the selection of tool variables to obtain the generalized moment estimator. Furthermore, the consistency of the obtained structural parameter estimators is proved by using the extreme value consistency theorem and the law of uniform large numbers (ULLN). To solve the problem of selecting spatial weight matrix, we can consider a crossvalidation method suitable for spatio-temporal data. Fortunately, a series of methods such as dimensionality reduction, feature selection, and model generalization has been provided by machine learning methods. The earliest crossvalidation method was called hold-out, which relied on only one partition of the data, and there was no crossover process, so it was also called the verification method [21]. Noting that the hold-out method relies on a partition of data and is easily affected by contingency factors, Geisser (2010) [22] proposed a crossvalidation method that includes the average of multiple hold-out estimates, realizing the transition from verification estimation to crossvalidation estimation. In order to reduce the combination number of data partition in crossvalidation, Shao (1993) [23] proposed the leave-P-out crossvalidation (LPOCV) in which the number of test samples in each data partition was the same. Especially in the special case when P=1, the method is evolved to leave-one-out crossvalidation (LOOCV). LOOCV is the simplest and most widely used crossvalidation in traditional analysis. Compared with the LPOCV considering all data partitioning, Geisser (2010) also proposed a crossvalidation based on only partial data partitioning, which is called RLT method. K-folded crossvalidation is proposed as an alternative to LOOCV which has a large computational overhead and relies on a basic partition of data divided into K-fold, each of which has a data capacity of N/K. In the case of limited samples, k-fold crossvalidation is the simplest and most widely used method of generalization error estimation. From the various crossvalidation methods that have appeared in the past, each method fully considers the randomness of the validation set to ensure the generalization ability of the test model. However, for the special panel data such as spatio-temporal data, there is usually an internal connection between spatial individuals, and the overall data also tends to have time trend. This problem is not taken into account by the previous crossvalidation methods, which may break the inherent regularity of spatio-temporal data. Based on the above considerations, this paper designs a kind of crossvalidation scheme suitable for spatio-temporal data. It is used to select stochastic frontier models, especially models with different weight matrices.Finally, the technology efficiency of China’s high-tech industry is analyzed by establishing a spatiotemporal stochastic frontier model.


## 2. Methodology

Previous studies on panel spatial stochastic frontier models mainly involved static panel spatial stochastic frontier models. Spatial lag effect is considered in the process of model building, but the influence of time lag effect is not included. If the time lag term and time-spatial lag term are added into the model, this kind of model can be called spatiotemporal stochastic frontier model. Obviously, the time-space double lag effect will produce stronger endogeneity, and new estimation methods should be considered to solve it.


### 2.1. Model Specification and Assumption

#### 2.1.1. Model Specification

The general form of the spatiotemporal stochastic frontier model can be stated in matrix form as(1)Yt=λ1WYt+λ2WYt‐1+γYt‐1+XtΒ+Εt,Et=ρMEt+εt,εt=vt−u,where Yt, Et, εt, and νt are N-dimensional vectors, whose components at time t=1,⋯,T are given by Yt=y1t,⋯,yNt′, Et=E1t,⋯,ENt′, εt=ε1t,⋯,εNt′, and vt=v1t,⋯,vNt′. The vector Yt consists of the outputs of the N production units, Et and εt are the composite error vectors corresponding to Yt, vt is the heterogeneous error vector, and u=u1,⋯,uN′ is the vector of time-invariant inefficiency terms. This kind of setting is appropriate when the time span is not large. As u is time invariant, it can be regarded as the individual effect, and thus, this paper primarily considers u as a fixed effect. Xt is an N×K-dimensional matrix consisting of the K exogenous input variables of the N production units at time t. W and M are N×N spatial weight matrices which are usually assumed to be different. If W=M, λ1 and ρ cannot be distinguished by means of the maximum likelihood method although they can be effectively distinguished by the GMM method [19]. The variables are stacked according to the section and time series in the following matrix form:
(2)Y=λ1IT⊗WY+λ2IT⊗WY‐1+γY‐1+XB+Ε,Ε=ρIT⊗MΕ+ε,ε=v−eT⊗INu,where ⊗ represents the Kronecker product of matrices, IT and IN are, respectively, the identity matrices of orders T and N, and eT is a T-dimensional column vector with all the entries equal to 1. The parameter vector of the model to be estimated is λ1,λ2,γ,Β,ρ,σv2,σu2, and its dimension is K+6, where Β is the parameters corresponding to the K-dimension exogenous explanatory variables Xt. Β and λ1,λ2,γ together constitute the structural parameters, and ρ,σv2,σu2 are the error term parameters.


#### 2.1.2. Model Assumption

The assumptions of the spatiotemporal stochastic frontier model are the following:Assumption 1.
The distribution of the error vectorv is given by v∼N0,σv2INT.Assumption 2.
The inefficiency termu is time invariant, with distribution u∼N+0,σu2INT.Assumption 3.
0<σv2<bv<∞, 0<σu2<bu<∞, and v and u have finite fourth moment.Assumption 4.
v and u are uncorrelated with X.Assumption 5.
The spatial weight matricesW and M satisfy wii=0 and mii=0; i=1,2,⋯,N. For arbitrary ∣λ1∣<1, ∣λ2∣<1, and ∣ρ∣<1, the matrices 1−λ1W, 1−λ2W, and 1−ρM are all nonsingular matrices. For each of the matrices W, M, 1−λ1W, 1−λ2W, 1−ρM, 1−λ1W−1, 1−λ2W−1, and 1−ρM−1, the row sums and column sums are all absolutely uniformly bounded.Assumption1 is a classic assumption of the spatial error autocorrelation model. By Assumption 2, the same individual inefficiency term remains constant at different times. When using GMM to estimate the structural parameters of the model, the distribution of error terms can be ignored; nevertheless, in order to improve the efficiency of computation, the half normal distribution for the inefficiency term is usually assumed. Assumption 3 ensures the boundedness of the variance of the error term in this model, which is an important condition for the consistency of the estimator. Assumption 4 is a classical assumption commonly used in traditional regression analysis methods, and the moment condition is set according to this assumption in the generalized moment estimation of this model. Assumption 5 is set according to the space weight matrix of this model and the properties of space station autoregressive coefficient and space-time autoregressive coefficient, which also ensures the consistency of parameter estimators.


### 2.2. Parameter Estimation

#### 2.2.1. Estimation strategy

There is an endogeneity problem in the model. For the spatial lag termλ1WYt in the model, there is
(3)COVWYt,Εt=COVWA−1λ2WYt−1+γYt−1+XtΒ+Εt,Εt=λ2COVWA−1WYt−1,Εt+γCOVWA−1Yt−1,Εt+COVWA−1Εt,Εt,where A=I−λ1W, and COVWA−1WYt−1,Εt=0,COVWA−1Yt−1,Εt=0, and COVWA−1XtΒ,Εt=0 can be obtained from the assumptions of the regression model, and COVWA−1Εt,Εt is quadratic. In spatial econometrics, W is usually not a zero matrix, and so, WA−1 is not a zero matrix. While taking into account the expected value of the compound error term cannot be 0, it can be considered that the quadratic form COVWA−1Εt,Εt is almost impossible to be equal to 0 (see Appendix A for proof). Therefore, in the dynamic panel spatial stochastic frontier models, there is an endogeneity problem which will lead to the inconsistency of traditional estimators. So, we considered GMM as a good way to solve the endogeneity problem.The parameter vector to be estimated in the model isλ1,λ2,γ,Β,ρ,σv2,σu2′, where λ1, λ2, γ, and Β are the structural parameters of the model, and ρ, σv2, and σu2 are the error distribution parameter of the model. The estimation of the model is completed in three steps:Step 1.
Using the GMM to estimate the structure parameterλ1,λ2,γ,Β in the model.Step 2.
Making a moment estimation of the parameterρ,σv2,σu2 that is included in the error term.Step 3.
Using the estimator obtained in Step2 to modify the result of Step 1.


#### 2.2.2. Estimation of Structural Parameterλ1,λ2,γ,Β

(1) Difference Model and Level Model. Anderson and Hsiao (1981) [24] proposed to use yi,t−2 as the instrumental variable of Δyi,t−1, and then, 2SLS estimation is carried out. This estimator is called “Anderson-Hsiao estimator.” According to the same logic, lag variables of higher order are also valid IV. Arellano and Bond (1991) [25] used all possible lag variables as IV (the number of IVs is more than the number of endogenous variables) to conduct GMM estimation. This GMM estimator is called Arellano-Bond estimator or difference GMM. The disadvantage of difference GMM is that the variable which does not change with time is eliminated, and its coefficient cannot be estimated. If the series yi,t has a strong persistence, that is, the first-order autoregressive coefficient is close to 1, then the correlation may be very weak and lead to the problem of weak instrumental variables. In order to solve the above two problems, Arellano and Bver (1995) [26] returned to the level equation and used Δyi,t−1,Δyi,t−2,⋯ as IV to estimate the GMM of the level equation, which was called “level GMM.” Blundell and Bond (1998) [27] combined difference GMM with level GMM and estimated the difference equation and level equation as one equation system for GMM, which was called “system GMM.” The advantage of system GMM is that it can improve the efficiency of estimation (small sample properties are better), and it can estimate the variable that does not change with time (the system GMM contains the level equation). In order to solve the endogenous problem of dynamic panel data model, Arellano and Bond (1991) [25], Arellano and Bover (1995) [26], and Blundell and Bond (1998) [27], respectively, considered from the perspective of difference model and level model, and different instrumental variables were selected.The corresponding difference model and level model of Equation (1) are simplified as
(4)ΔYt=ΔZtθ+ΔΕt,(5)Yt=Ztθ+Εt.(4) and (5) can also be collectively called spatial system model, where Equation (4) is the difference model, and Equation (5) is the level model, Zt=WYt,WYt−1,Yt−1,Xt′ is the vector composed of all explanatory variables, and θ=λ1,λ2,γ,B′ is the vector composed of structural parameters. The expansion of Equation (5) is Equation (1); the expansion of Equation (4) can be expressed as follows:
(6)ΔYt=λ1WΔYt+λ2WΔYt‐1+γΔYt‐1+ΔXtΒ+ΔΕt,Et=ρMΔEt+Δεt,Δεt=Δvt.(2) Moment Condition and Instrumental Variable. Since ΔXt is a strictly exogenous variable, it is not related to the compound error term ΔΕt, nor is it related to Εt. The moment conditions for identifying B in the difference model and the level model are as follows:
(7)EΔXt′ΔΕt=0,t=3,⋯,T,EΔXt′Εt=0,t=3,⋯,T.The moment condition structure for identifyingλ2 and γ in the two models is as follows: since the spatial lag term and time lag term of the dependent variable ΔYt are both endogenous variables, therefore, it is necessary to find a set of instrumental variables that is related to time lag and space lag and exogenous explanatory variables, but not related to the difference error term ΔΕtt=3,⋯,T. Arellano and Bond (1991) [25] uses all possible level lag variables yt−2,⋯,y1 of Yt as instrumental variables for the time-lag first-order difference term (ΔYt−1) of the dependent variable. These instrumental variables are related to (ΔYt−1), but not to ΔΕt. The moment conditions corresponding to the difference model and the level model are as follows:
(8)EYt−s′ΔΕt=0,t=3,⋯,T;s=2,⋯,t−1,(9)EΔYt−s′Εt=0,t=3,⋯,T;s=2,⋯,t−1.The moment conditions for identifyingλ1 in the two models are as follows.Construct a spatial lag itemWYt as follows; Jacobs et al. (2009) [19] provided a method of finding instrumental variables, that is time lag terms of spatial lag dependent variables, who also proved that the moment condition obtained by this method was as valid as Equation (8). So, corresponding to the difference model and the level model, the following moment conditions can be listed:
(10)EWlYt−s′ΔΕt=0,t=3,⋯,T;s=2,⋯,t−1;l=1,⋯,L,EWlΔYt−s′Εt=0,t=3,⋯,T;s=1,⋯,t−2;l=1,⋯L,where l is the exponential of matrix W and the integer L is the maximum order of spatial lag that can be used as the instrumental variable.In addition, based on the method provided by Kelejian and Robinson (1993) [28], formula (1) shows that WYt depends on WXt, so the instrumental variable WΔXt can be selected by the first-order difference method for WΔYt. Since ΔXt is a strictly exogenous variable, it is not related to the compound error term ΔEt, so corresponding to the difference model and the level model, the instrumental variables satisfy the following moment conditions:
(11)EWlΔXt′ΔΕt=0,t=3,⋯,T,EWlΔXt′Εt=0,t=3,⋯,T.(3) GMM Estimation. When we estimate the parameters of the spatio-temporal stochastic frontier model, we use the system GMM method similar to the general dynamic panel model to construct the spatial system GMM estimation. Unlike the system GMM, the IVs of the spatial system GMM are composed of time lag variable and spatial lag variable.For each period oft, the moment condition of J≥K+2 can be given. The moment conditions corresponding to the difference model and the level model can be abbreviated as
(12)EHN,ABt′ΔΕt=0,EHN,Lt′Εt=0.The matricesHN,ABt and HN,Lt are expressed as follows
(13)Yt−s=y100000⋯⋯00⋯00y2y1000⋯⋯00⋯0000y3y2y1⋯⋯00⋯0⋮⋮⋮⋮⋮⋮⋯⋯00⋯0000000⋯⋯yt−2yt−3⋯y1,Xt=x3x4⋮xt,ΔYt−s=Δy100000⋯⋯00⋯00Δy2Δy1000⋯⋯00⋯0000Δy3Δy2Δy1⋯⋯00⋯0⋮⋮⋮⋮⋮⋮⋯⋯00⋯0000000⋯⋯Δyt−2Δyt−3⋯Δy1,ΔXt=Δx3Δx4⋮Δxt.Then,HN,ABt=Yt−s,WYt−s,W2Yt−s,⋯,WLYt−s,WΔXt,ΔXt and
(14)HN,Lt=ΔYt−s,WΔYt−s,W2ΔYt−s,⋯,WLΔYt−s,WXt,Xtt=3,⋯,T.That is,HN,ABt and HN,Lt are matrices with instrumental variables as column vectors, and the subscript N means that the matrix depends on the unit number of individuals. Let HN,AB and HN,L be block diagonal matrices composed of block HN,ABt and HN,Lt, respectively. In order to define the GMM (Spatial Blundell Bond, SBB) estimator of the spatial dynamic panel stochastic frontier model, the difference variables and level variables are combined to define the matrix as follows:
(15)Y∗N=ΔY′N,Y′N,Z∗N=ΔZ′N,Z′N,Ε∗N=ΔΕ′N,Ε′N.The instrument variable matrix is(16)HN,BB=diagHN,AB,HN,L=HN,AB00HN,L,where HN,AB is the instrumental matrix of spatial difference GMM estimation, and HN,L is the instrumental matrix of spatial level GMM estimation. The weight matrix is
(17)GN,BB=diagGN,AB,IT−2⊗IN=GN,AB00IT−2⊗IN.This diagonal of the matrix is composed of the weight matrix defined in the process of spatial difference GMM estimation and an identity matrix, whereGN,AB=IN⊗G is NT−2×NT−2 weight matrix which elements are
(18)Gij=2,i=j,−1,i=j+1,−1,j=i+1,0,others.This weight matrix is proposed by Arellano and Bond (1991) [25] which is further define the weight matrix:
(19)AN,BB=H′N,BBGN,BBHN,BBN−1.Through the above process, combining the spatial difference equation with the spatial level equation, we get the spatial system GMM estimation process. Get the objective function of generalized moment estimation for spatial system as follows:(20)1NH′N,BBΔΕ∗N′AN,BBH′N,BBΔΕ∗N=1NH′N,BBΔY∗N−ΔZ∗NθN′AN,BBH′N,BBΔY∗N−ΔZ∗NθN.The one-stage SBB estimator ofθ can be obtained by minimizing Equation (20):
(21)θ^SBB=Z∗′NH′N,BBA∗NH′N,BBZ∗NZ∗′NH′N,BBA∗NH′N,BBY∗N.Equation (21) can also be called the spatial system GMM estimator.(4) Improvement of Instrumental Variable Matrix. The instrumental variable matrix constructed in accordance with the above method has a high dimension and grows exponentially as the values of T and L increase. In order to reduce the dimension of the instrumental variable matrix and avoid overfitting the instrumental variable, we can simplify it by using the “condensing instrumental variable matrix” proposed by Beck and Levine (2004) [29].We still sets=2 and L=1 in the GMM instrumental variable matrix of the space system, and the corresponding condensed instrumental variable matrix is
(22)H′N,BB=H1AB⋯HNABH1L⋯HNL,where HiL (i=1,⋯,N) is the instrumental variable quantum matrix of the level model corresponding to the i individual.
(23)HiL=Δy100⋯WΔy100⋯WΔX3ΔX3Δy2Δy10⋯WΔy2WΔy10⋯WΔX4ΔX4Δy3Δy2Δy1⋯WΔy3WΔy2WΔy1⋯WΔX5ΔX5⋮⋮⋮⋱⋮⋮⋮⋱⋮⋮.


### 2.3. Estimation of Error Term Parameters

#### 2.3.1. Estimation ofρ,σv2

After the estimator of the structural parameterθN=λ1,λ2,γ,Β′ in model (1) is obtained in the first stage, the model residual Ε^t=Yt−Ztθ^N can be further obtained, in which the Zt=WYt,WYt−1,Yt−1,Xt is the vector set composed of all explanatory variables in model (1). Consistent GMM estimation can be obtained by using residual Ε^t and modifying the moment condition proposed by Kapoor et al. (2007) [18]. The specific process is as follows.According to the assumptions of the model, the individual effect of the model is the inefficiency term. According to the covariance structure of the compound error term, it can be known that(24)Covεε′=σu2IN⊗JT+σv2INT.Introduce transformation matrix:(25)Q0=IT−JTT⊗IN,Q1=JTT⊗IN,where IT and IN are the identity matrix of order T and order N, JT=eTeT′ is the matrix of order T×T, and the elements of that are all 1. Properties of transformation matrix Q0: (i) trQ0=NT−1; (ii) Q0eT⊗IN=0; (iii) Q0IT⊗RN=IT⊗RNQ0 (where RN is any N×N matrix). From the properties (ii), we can further deduce the special properties (iv) Q0ε=Q0v of matrix Q0 in this paper, where ε and ν are the corresponding error terms in model (1).Let(26)Ε¯=IT⊗MΕ;Ε¯¯=IT⊗MΕ¯;ε¯=IT⊗Mε.Then,(27)ε=Ε−ρΕ¯,ε¯=Ε¯−ρΕ¯¯.Based on the above transformation and referring to the first three of the six moment conditions given by Kapoor et al. (2007) [18] and related properties, the following three moment conditions are given in this paper:
(28)Eε′Q0ε=σv2NT−1,Eε¯′Q0ε¯=σv2NT−1trM′M,Eε¯′Q0ε=0.To further integrate the above moment conditions, we can get that(29)E1NT−1ε′Q0ε1NT−1ε¯′Q0ε¯1NT−1ε¯′Q0ε=σv21Nσv2trM′M0.Substitute Equations (26) and (27) into Equation (29) to obtain that
(30)1NE2T−1Ε¯′Q0E−1T−1Ε¯′Q0Ε¯12T−1Ε¯′Q0Ε¯¯−1T−1Ε¯¯′Q0Ε¯¯trM′M1T−1Ε¯′Q0Ε¯¯+Ε¯′Q0Ε¯−1T−1Ε¯′Q0Ε¯¯0ρρ2σv2=1N1T−1E′Q0E1T−1Ε¯′Q0Ε¯1T−1E′Q0Ε¯.The residualΕ^t=Yt−Ztθ^N estimated in the first stage is substituted into Ε¯ and Ε¯¯ in Equation (30) to obtain the sample moment equation. In the sample moment equation, the estimated value ρ^,σ^v2 of ρ and σv2 can be solved by the following objective function:
(31)ρ^,σ^v2=argminGNρ,ρ2,σv2−gN′GNρ,ρ2,σv2−gN,ρ∈−a,a,σv2∈0,bv,where
(32)GN=1N2T−1Ε∧¯′Q0Ε^−1T−1Ε∧¯′Q0Ε^¯12T−1Ε∧¯′Q0Ε^¯¯−1T−1Ε∧¯¯′Q0Ε^¯¯trM′M1T−1Ε∧¯′Q0Ε^¯¯+Ε∧¯′Q0Ε^¯−1T−1Ε∧¯′Q0Ε^¯¯0,gN=1N1T−1Ε∧′Q0Ε^1T−1Ε∧¯′Q0Ε^¯1T−1Ε∧′Q0Ε^¯.


#### 2.3.2. Estimation ofσu2

The fourth moment condition given by Kapoor et al. (2007) [18] is
(33)Eε′Q1ε=NTσu2+Nσv2=Nσ12,where σ12=Tσu2+σv2. However, considering the characteristics of the stochastic frontier model, the compound error term ε obeys the asymmetric distribution of the expected nonzero, so the moment condition (33) cannot be directly applied, and the following formula can be proved:
(34)COVε′Q1ε=NTσu2+Nσv2=Nσ12.The estimator of parameterσ12 can be obtained by combining Equation (34) with Equation (31):
(35)σ^12=1NCOVΕ∧−ρ∧Ε∧¯′Q1Ε^−ρ^Ε^¯=1T−1COVΕ∧′Q1Ε^−2T−1COVΕ∧¯′Q1Ε^·ρ^−1T−1COVΕ∧¯′Q1Ε^¯ρ∧2.Substituteσ^v2 and σ^12 into (33) to obtain the moment estimator σ^u2 of σu2.


### 2.4. Spatial Correction of Estimators

Although it can be proved that the estimator (21) is a consistent estimator, it can also be proved that the consistency of the GMM estimator can be guaranteed even if the model has spatial error autocorrelation. However, the estimator (21) cannot solve the spatial dependence of the error term, and the variance of the estimator is relatively large. After obtaining the consistent estimator of ρ by Equation (31), the consistent estimator can be obtained by a correcting transformation. According to the spatial correction method given by Jacobs et al. (2009) [19], the estimator obtained in the first step was corrected.The estimatorρ^ obtained from Equation (31) is used to construct matrix I−ρ^M, and left the difference GMM and the explained variables and the instrumental variables matrix of the system GMM estimation, if
(36)Y~=I−ρ^MY,ΔY~N=I−ρ^MΔYN,Y~−1=I−ρ^MY−1,W~Y=I−ρ^MWY,X~=I−ρ^MX.The corresponding explanatory variable setZt=WYt,WYt−1,Yt−1,Xt is corrected as
(37)Z~=I−ρ^MZ∗N,ΔZ~=I−ρ^MΔZN.The instrumental variable matrix and weight matrix corresponding to the GMM estimation of the spatial system are corrected as follows:(38)H~SBB=I−ρ^MHN,BB,A~SBB=H~′SBBH~SBB−1.Then, the corrected system GMM estimator is(39)θ~SBB=Z~′SBBH~SBBA~SBBH~′SBBZ~SBBZ~′SBBH~SBBA~SBBH~′SBBY~SBB.


## 2.1. Model Specification and Assumption

### 2.1.1. Model Specification

The general form of the spatiotemporal stochastic frontier model can be stated in matrix form as(1)Yt=λ1WYt+λ2WYt‐1+γYt‐1+XtΒ+Εt,Et=ρMEt+εt,εt=vt−u,where Yt, Et, εt, and νt are N-dimensional vectors, whose components at time t=1,⋯,T are given by Yt=y1t,⋯,yNt′, Et=E1t,⋯,ENt′, εt=ε1t,⋯,εNt′, and vt=v1t,⋯,vNt′. The vector Yt consists of the outputs of the N production units, Et and εt are the composite error vectors corresponding to Yt, vt is the heterogeneous error vector, and u=u1,⋯,uN′ is the vector of time-invariant inefficiency terms. This kind of setting is appropriate when the time span is not large. As u is time invariant, it can be regarded as the individual effect, and thus, this paper primarily considers u as a fixed effect. Xt is an N×K-dimensional matrix consisting of the K exogenous input variables of the N production units at time t. W and M are N×N spatial weight matrices which are usually assumed to be different. If W=M, λ1 and ρ cannot be distinguished by means of the maximum likelihood method although they can be effectively distinguished by the GMM method [19]. The variables are stacked according to the section and time series in the following matrix form:
(2)Y=λ1IT⊗WY+λ2IT⊗WY‐1+γY‐1+XB+Ε,Ε=ρIT⊗MΕ+ε,ε=v−eT⊗INu,where ⊗ represents the Kronecker product of matrices, IT and IN are, respectively, the identity matrices of orders T and N, and eT is a T-dimensional column vector with all the entries equal to 1. The parameter vector of the model to be estimated is λ1,λ2,γ,Β,ρ,σv2,σu2, and its dimension is K+6, where Β is the parameters corresponding to the K-dimension exogenous explanatory variables Xt. Β and λ1,λ2,γ together constitute the structural parameters, and ρ,σv2,σu2 are the error term parameters.


### 2.1.2. Model Assumption

The assumptions of the spatiotemporal stochastic frontier model are the following:Assumption 1.
The distribution of the error vectorv is given by v∼N0,σv2INT.Assumption 2.
The inefficiency termu is time invariant, with distribution u∼N+0,σu2INT.Assumption 3.
0<σv2<bv<∞, 0<σu2<bu<∞, and v and u have finite fourth moment.Assumption 4.
v and u are uncorrelated with X.Assumption 5.
The spatial weight matricesW and M satisfy wii=0 and mii=0; i=1,2,⋯,N. For arbitrary ∣λ1∣<1, ∣λ2∣<1, and ∣ρ∣<1, the matrices 1−λ1W, 1−λ2W, and 1−ρM are all nonsingular matrices. For each of the matrices W, M, 1−λ1W, 1−λ2W, 1−ρM, 1−λ1W−1, 1−λ2W−1, and 1−ρM−1, the row sums and column sums are all absolutely uniformly bounded.Assumption1 is a classic assumption of the spatial error autocorrelation model. By Assumption 2, the same individual inefficiency term remains constant at different times. When using GMM to estimate the structural parameters of the model, the distribution of error terms can be ignored; nevertheless, in order to improve the efficiency of computation, the half normal distribution for the inefficiency term is usually assumed. Assumption 3 ensures the boundedness of the variance of the error term in this model, which is an important condition for the consistency of the estimator. Assumption 4 is a classical assumption commonly used in traditional regression analysis methods, and the moment condition is set according to this assumption in the generalized moment estimation of this model. Assumption 5 is set according to the space weight matrix of this model and the properties of space station autoregressive coefficient and space-time autoregressive coefficient, which also ensures the consistency of parameter estimators.


## 2.1.1. Model Specification

The general form of the spatiotemporal stochastic frontier model can be stated in matrix form as(1)Yt=λ1WYt+λ2WYt‐1+γYt‐1+XtΒ+Εt,Et=ρMEt+εt,εt=vt−u,where Yt, Et, εt, and νt are N-dimensional vectors, whose components at time t=1,⋯,T are given by Yt=y1t,⋯,yNt′, Et=E1t,⋯,ENt′, εt=ε1t,⋯,εNt′, and vt=v1t,⋯,vNt′. The vector Yt consists of the outputs of the N production units, Et and εt are the composite error vectors corresponding to Yt, vt is the heterogeneous error vector, and u=u1,⋯,uN′ is the vector of time-invariant inefficiency terms. This kind of setting is appropriate when the time span is not large. As u is time invariant, it can be regarded as the individual effect, and thus, this paper primarily considers u as a fixed effect. Xt is an N×K-dimensional matrix consisting of the K exogenous input variables of the N production units at time t. W and M are N×N spatial weight matrices which are usually assumed to be different. If W=M, λ1 and ρ cannot be distinguished by means of the maximum likelihood method although they can be effectively distinguished by the GMM method [19]. The variables are stacked according to the section and time series in the following matrix form:
(2)Y=λ1IT⊗WY+λ2IT⊗WY‐1+γY‐1+XB+Ε,Ε=ρIT⊗MΕ+ε,ε=v−eT⊗INu,where ⊗ represents the Kronecker product of matrices, IT and IN are, respectively, the identity matrices of orders T and N, and eT is a T-dimensional column vector with all the entries equal to 1. The parameter vector of the model to be estimated is λ1,λ2,γ,Β,ρ,σv2,σu2, and its dimension is K+6, where Β is the parameters corresponding to the K-dimension exogenous explanatory variables Xt. Β and λ1,λ2,γ together constitute the structural parameters, and ρ,σv2,σu2 are the error term parameters.


## 2.1.2. Model Assumption

The assumptions of the spatiotemporal stochastic frontier model are the following:Assumption 1.
The distribution of the error vectorv is given by v∼N0,σv2INT.Assumption 2.
The inefficiency termu is time invariant, with distribution u∼N+0,σu2INT.Assumption 3.
0<σv2<bv<∞, 0<σu2<bu<∞, and v and u have finite fourth moment.Assumption 4.
v and u are uncorrelated with X.Assumption 5.
The spatial weight matricesW and M satisfy wii=0 and mii=0; i=1,2,⋯,N. For arbitrary ∣λ1∣<1, ∣λ2∣<1, and ∣ρ∣<1, the matrices 1−λ1W, 1−λ2W, and 1−ρM are all nonsingular matrices. For each of the matrices W, M, 1−λ1W, 1−λ2W, 1−ρM, 1−λ1W−1, 1−λ2W−1, and 1−ρM−1, the row sums and column sums are all absolutely uniformly bounded.Assumption1 is a classic assumption of the spatial error autocorrelation model. By Assumption 2, the same individual inefficiency term remains constant at different times. When using GMM to estimate the structural parameters of the model, the distribution of error terms can be ignored; nevertheless, in order to improve the efficiency of computation, the half normal distribution for the inefficiency term is usually assumed. Assumption 3 ensures the boundedness of the variance of the error term in this model, which is an important condition for the consistency of the estimator. Assumption 4 is a classical assumption commonly used in traditional regression analysis methods, and the moment condition is set according to this assumption in the generalized moment estimation of this model. Assumption 5 is set according to the space weight matrix of this model and the properties of space station autoregressive coefficient and space-time autoregressive coefficient, which also ensures the consistency of parameter estimators.


## 2.2. Parameter Estimation

### 2.2.1. Estimation strategy

There is an endogeneity problem in the model. For the spatial lag termλ1WYt in the model, there is
(3)COVWYt,Εt=COVWA−1λ2WYt−1+γYt−1+XtΒ+Εt,Εt=λ2COVWA−1WYt−1,Εt+γCOVWA−1Yt−1,Εt+COVWA−1Εt,Εt,where A=I−λ1W, and COVWA−1WYt−1,Εt=0,COVWA−1Yt−1,Εt=0, and COVWA−1XtΒ,Εt=0 can be obtained from the assumptions of the regression model, and COVWA−1Εt,Εt is quadratic. In spatial econometrics, W is usually not a zero matrix, and so, WA−1 is not a zero matrix. While taking into account the expected value of the compound error term cannot be 0, it can be considered that the quadratic form COVWA−1Εt,Εt is almost impossible to be equal to 0 (see Appendix A for proof). Therefore, in the dynamic panel spatial stochastic frontier models, there is an endogeneity problem which will lead to the inconsistency of traditional estimators. So, we considered GMM as a good way to solve the endogeneity problem.The parameter vector to be estimated in the model isλ1,λ2,γ,Β,ρ,σv2,σu2′, where λ1, λ2, γ, and Β are the structural parameters of the model, and ρ, σv2, and σu2 are the error distribution parameter of the model. The estimation of the model is completed in three steps:Step 1.
Using the GMM to estimate the structure parameterλ1,λ2,γ,Β in the model.Step 2.
Making a moment estimation of the parameterρ,σv2,σu2 that is included in the error term.Step 3.
Using the estimator obtained in Step2 to modify the result of Step 1.


### 2.2.2. Estimation of Structural Parameterλ1,λ2,γ,Β

(1) Difference Model and Level Model. Anderson and Hsiao (1981) [24] proposed to use yi,t−2 as the instrumental variable of Δyi,t−1, and then, 2SLS estimation is carried out. This estimator is called “Anderson-Hsiao estimator.” According to the same logic, lag variables of higher order are also valid IV. Arellano and Bond (1991) [25] used all possible lag variables as IV (the number of IVs is more than the number of endogenous variables) to conduct GMM estimation. This GMM estimator is called Arellano-Bond estimator or difference GMM. The disadvantage of difference GMM is that the variable which does not change with time is eliminated, and its coefficient cannot be estimated. If the series yi,t has a strong persistence, that is, the first-order autoregressive coefficient is close to 1, then the correlation may be very weak and lead to the problem of weak instrumental variables. In order to solve the above two problems, Arellano and Bver (1995) [26] returned to the level equation and used Δyi,t−1,Δyi,t−2,⋯ as IV to estimate the GMM of the level equation, which was called “level GMM.” Blundell and Bond (1998) [27] combined difference GMM with level GMM and estimated the difference equation and level equation as one equation system for GMM, which was called “system GMM.” The advantage of system GMM is that it can improve the efficiency of estimation (small sample properties are better), and it can estimate the variable that does not change with time (the system GMM contains the level equation). In order to solve the endogenous problem of dynamic panel data model, Arellano and Bond (1991) [25], Arellano and Bover (1995) [26], and Blundell and Bond (1998) [27], respectively, considered from the perspective of difference model and level model, and different instrumental variables were selected.The corresponding difference model and level model of Equation (1) are simplified as
(4)ΔYt=ΔZtθ+ΔΕt,(5)Yt=Ztθ+Εt.(4) and (5) can also be collectively called spatial system model, where Equation (4) is the difference model, and Equation (5) is the level model, Zt=WYt,WYt−1,Yt−1,Xt′ is the vector composed of all explanatory variables, and θ=λ1,λ2,γ,B′ is the vector composed of structural parameters. The expansion of Equation (5) is Equation (1); the expansion of Equation (4) can be expressed as follows:
(6)ΔYt=λ1WΔYt+λ2WΔYt‐1+γΔYt‐1+ΔXtΒ+ΔΕt,Et=ρMΔEt+Δεt,Δεt=Δvt.(2) Moment Condition and Instrumental Variable. Since ΔXt is a strictly exogenous variable, it is not related to the compound error term ΔΕt, nor is it related to Εt. The moment conditions for identifying B in the difference model and the level model are as follows:
(7)EΔXt′ΔΕt=0,t=3,⋯,T,EΔXt′Εt=0,t=3,⋯,T.The moment condition structure for identifyingλ2 and γ in the two models is as follows: since the spatial lag term and time lag term of the dependent variable ΔYt are both endogenous variables, therefore, it is necessary to find a set of instrumental variables that is related to time lag and space lag and exogenous explanatory variables, but not related to the difference error term ΔΕtt=3,⋯,T. Arellano and Bond (1991) [25] uses all possible level lag variables yt−2,⋯,y1 of Yt as instrumental variables for the time-lag first-order difference term (ΔYt−1) of the dependent variable. These instrumental variables are related to (ΔYt−1), but not to ΔΕt. The moment conditions corresponding to the difference model and the level model are as follows:
(8)EYt−s′ΔΕt=0,t=3,⋯,T;s=2,⋯,t−1,(9)EΔYt−s′Εt=0,t=3,⋯,T;s=2,⋯,t−1.The moment conditions for identifyingλ1 in the two models are as follows.Construct a spatial lag itemWYt as follows; Jacobs et al. (2009) [19] provided a method of finding instrumental variables, that is time lag terms of spatial lag dependent variables, who also proved that the moment condition obtained by this method was as valid as Equation (8). So, corresponding to the difference model and the level model, the following moment conditions can be listed:
(10)EWlYt−s′ΔΕt=0,t=3,⋯,T;s=2,⋯,t−1;l=1,⋯,L,EWlΔYt−s′Εt=0,t=3,⋯,T;s=1,⋯,t−2;l=1,⋯L,where l is the exponential of matrix W and the integer L is the maximum order of spatial lag that can be used as the instrumental variable.In addition, based on the method provided by Kelejian and Robinson (1993) [28], formula (1) shows that WYt depends on WXt, so the instrumental variable WΔXt can be selected by the first-order difference method for WΔYt. Since ΔXt is a strictly exogenous variable, it is not related to the compound error term ΔEt, so corresponding to the difference model and the level model, the instrumental variables satisfy the following moment conditions:
(11)EWlΔXt′ΔΕt=0,t=3,⋯,T,EWlΔXt′Εt=0,t=3,⋯,T.(3) GMM Estimation. When we estimate the parameters of the spatio-temporal stochastic frontier model, we use the system GMM method similar to the general dynamic panel model to construct the spatial system GMM estimation. Unlike the system GMM, the IVs of the spatial system GMM are composed of time lag variable and spatial lag variable.For each period oft, the moment condition of J≥K+2 can be given. The moment conditions corresponding to the difference model and the level model can be abbreviated as
(12)EHN,ABt′ΔΕt=0,EHN,Lt′Εt=0.The matricesHN,ABt and HN,Lt are expressed as follows
(13)Yt−s=y100000⋯⋯00⋯00y2y1000⋯⋯00⋯0000y3y2y1⋯⋯00⋯0⋮⋮⋮⋮⋮⋮⋯⋯00⋯0000000⋯⋯yt−2yt−3⋯y1,Xt=x3x4⋮xt,ΔYt−s=Δy100000⋯⋯00⋯00Δy2Δy1000⋯⋯00⋯0000Δy3Δy2Δy1⋯⋯00⋯0⋮⋮⋮⋮⋮⋮⋯⋯00⋯0000000⋯⋯Δyt−2Δyt−3⋯Δy1,ΔXt=Δx3Δx4⋮Δxt.Then,HN,ABt=Yt−s,WYt−s,W2Yt−s,⋯,WLYt−s,WΔXt,ΔXt and
(14)HN,Lt=ΔYt−s,WΔYt−s,W2ΔYt−s,⋯,WLΔYt−s,WXt,Xtt=3,⋯,T.That is,HN,ABt and HN,Lt are matrices with instrumental variables as column vectors, and the subscript N means that the matrix depends on the unit number of individuals. Let HN,AB and HN,L be block diagonal matrices composed of block HN,ABt and HN,Lt, respectively. In order to define the GMM (Spatial Blundell Bond, SBB) estimator of the spatial dynamic panel stochastic frontier model, the difference variables and level variables are combined to define the matrix as follows:
(15)Y∗N=ΔY′N,Y′N,Z∗N=ΔZ′N,Z′N,Ε∗N=ΔΕ′N,Ε′N.The instrument variable matrix is(16)HN,BB=diagHN,AB,HN,L=HN,AB00HN,L,where HN,AB is the instrumental matrix of spatial difference GMM estimation, and HN,L is the instrumental matrix of spatial level GMM estimation. The weight matrix is
(17)GN,BB=diagGN,AB,IT−2⊗IN=GN,AB00IT−2⊗IN.This diagonal of the matrix is composed of the weight matrix defined in the process of spatial difference GMM estimation and an identity matrix, whereGN,AB=IN⊗G is NT−2×NT−2 weight matrix which elements are
(18)Gij=2,i=j,−1,i=j+1,−1,j=i+1,0,others.This weight matrix is proposed by Arellano and Bond (1991) [25] which is further define the weight matrix:
(19)AN,BB=H′N,BBGN,BBHN,BBN−1.Through the above process, combining the spatial difference equation with the spatial level equation, we get the spatial system GMM estimation process. Get the objective function of generalized moment estimation for spatial system as follows:(20)1NH′N,BBΔΕ∗N′AN,BBH′N,BBΔΕ∗N=1NH′N,BBΔY∗N−ΔZ∗NθN′AN,BBH′N,BBΔY∗N−ΔZ∗NθN.The one-stage SBB estimator ofθ can be obtained by minimizing Equation (20):
(21)θ^SBB=Z∗′NH′N,BBA∗NH′N,BBZ∗NZ∗′NH′N,BBA∗NH′N,BBY∗N.Equation (21) can also be called the spatial system GMM estimator.(4) Improvement of Instrumental Variable Matrix. The instrumental variable matrix constructed in accordance with the above method has a high dimension and grows exponentially as the values of T and L increase. In order to reduce the dimension of the instrumental variable matrix and avoid overfitting the instrumental variable, we can simplify it by using the “condensing instrumental variable matrix” proposed by Beck and Levine (2004) [29].We still sets=2 and L=1 in the GMM instrumental variable matrix of the space system, and the corresponding condensed instrumental variable matrix is
(22)H′N,BB=H1AB⋯HNABH1L⋯HNL,where HiL (i=1,⋯,N) is the instrumental variable quantum matrix of the level model corresponding to the i individual.
(23)HiL=Δy100⋯WΔy100⋯WΔX3ΔX3Δy2Δy10⋯WΔy2WΔy10⋯WΔX4ΔX4Δy3Δy2Δy1⋯WΔy3WΔy2WΔy1⋯WΔX5ΔX5⋮⋮⋮⋱⋮⋮⋮⋱⋮⋮.


## 2.2.1. Estimation strategy

There is an endogeneity problem in the model. For the spatial lag termλ1WYt in the model, there is
(3)COVWYt,Εt=COVWA−1λ2WYt−1+γYt−1+XtΒ+Εt,Εt=λ2COVWA−1WYt−1,Εt+γCOVWA−1Yt−1,Εt+COVWA−1Εt,Εt,where A=I−λ1W, and COVWA−1WYt−1,Εt=0,COVWA−1Yt−1,Εt=0, and COVWA−1XtΒ,Εt=0 can be obtained from the assumptions of the regression model, and COVWA−1Εt,Εt is quadratic. In spatial econometrics, W is usually not a zero matrix, and so, WA−1 is not a zero matrix. While taking into account the expected value of the compound error term cannot be 0, it can be considered that the quadratic form COVWA−1Εt,Εt is almost impossible to be equal to 0 (see Appendix A for proof). Therefore, in the dynamic panel spatial stochastic frontier models, there is an endogeneity problem which will lead to the inconsistency of traditional estimators. So, we considered GMM as a good way to solve the endogeneity problem.The parameter vector to be estimated in the model isλ1,λ2,γ,Β,ρ,σv2,σu2′, where λ1, λ2, γ, and Β are the structural parameters of the model, and ρ, σv2, and σu2 are the error distribution parameter of the model. The estimation of the model is completed in three steps:Step 1.
Using the GMM to estimate the structure parameterλ1,λ2,γ,Β in the model.Step 2.
Making a moment estimation of the parameterρ,σv2,σu2 that is included in the error term.Step 3.
Using the estimator obtained in Step2 to modify the result of Step 1.


## 2.2.2. Estimation of Structural Parameterλ1,λ2,γ,Β

(1) Difference Model and Level Model. Anderson and Hsiao (1981) [24] proposed to use yi,t−2 as the instrumental variable of Δyi,t−1, and then, 2SLS estimation is carried out. This estimator is called “Anderson-Hsiao estimator.” According to the same logic, lag variables of higher order are also valid IV. Arellano and Bond (1991) [25] used all possible lag variables as IV (the number of IVs is more than the number of endogenous variables) to conduct GMM estimation. This GMM estimator is called Arellano-Bond estimator or difference GMM. The disadvantage of difference GMM is that the variable which does not change with time is eliminated, and its coefficient cannot be estimated. If the series yi,t has a strong persistence, that is, the first-order autoregressive coefficient is close to 1, then the correlation may be very weak and lead to the problem of weak instrumental variables. In order to solve the above two problems, Arellano and Bver (1995) [26] returned to the level equation and used Δyi,t−1,Δyi,t−2,⋯ as IV to estimate the GMM of the level equation, which was called “level GMM.” Blundell and Bond (1998) [27] combined difference GMM with level GMM and estimated the difference equation and level equation as one equation system for GMM, which was called “system GMM.” The advantage of system GMM is that it can improve the efficiency of estimation (small sample properties are better), and it can estimate the variable that does not change with time (the system GMM contains the level equation). In order to solve the endogenous problem of dynamic panel data model, Arellano and Bond (1991) [25], Arellano and Bover (1995) [26], and Blundell and Bond (1998) [27], respectively, considered from the perspective of difference model and level model, and different instrumental variables were selected.The corresponding difference model and level model of Equation (1) are simplified as
(4)ΔYt=ΔZtθ+ΔΕt,(5)Yt=Ztθ+Εt.(4) and (5) can also be collectively called spatial system model, where Equation (4) is the difference model, and Equation (5) is the level model, Zt=WYt,WYt−1,Yt−1,Xt′ is the vector composed of all explanatory variables, and θ=λ1,λ2,γ,B′ is the vector composed of structural parameters. The expansion of Equation (5) is Equation (1); the expansion of Equation (4) can be expressed as follows:
(6)ΔYt=λ1WΔYt+λ2WΔYt‐1+γΔYt‐1+ΔXtΒ+ΔΕt,Et=ρMΔEt+Δεt,Δεt=Δvt.(2) Moment Condition and Instrumental Variable. Since ΔXt is a strictly exogenous variable, it is not related to the compound error term ΔΕt, nor is it related to Εt. The moment conditions for identifying B in the difference model and the level model are as follows:
(7)EΔXt′ΔΕt=0,t=3,⋯,T,EΔXt′Εt=0,t=3,⋯,T.The moment condition structure for identifyingλ2 and γ in the two models is as follows: since the spatial lag term and time lag term of the dependent variable ΔYt are both endogenous variables, therefore, it is necessary to find a set of instrumental variables that is related to time lag and space lag and exogenous explanatory variables, but not related to the difference error term ΔΕtt=3,⋯,T. Arellano and Bond (1991) [25] uses all possible level lag variables yt−2,⋯,y1 of Yt as instrumental variables for the time-lag first-order difference term (ΔYt−1) of the dependent variable. These instrumental variables are related to (ΔYt−1), but not to ΔΕt. The moment conditions corresponding to the difference model and the level model are as follows:
(8)EYt−s′ΔΕt=0,t=3,⋯,T;s=2,⋯,t−1,(9)EΔYt−s′Εt=0,t=3,⋯,T;s=2,⋯,t−1.The moment conditions for identifyingλ1 in the two models are as follows.Construct a spatial lag itemWYt as follows; Jacobs et al. (2009) [19] provided a method of finding instrumental variables, that is time lag terms of spatial lag dependent variables, who also proved that the moment condition obtained by this method was as valid as Equation (8). So, corresponding to the difference model and the level model, the following moment conditions can be listed:
(10)EWlYt−s′ΔΕt=0,t=3,⋯,T;s=2,⋯,t−1;l=1,⋯,L,EWlΔYt−s′Εt=0,t=3,⋯,T;s=1,⋯,t−2;l=1,⋯L,where l is the exponential of matrix W and the integer L is the maximum order of spatial lag that can be used as the instrumental variable.In addition, based on the method provided by Kelejian and Robinson (1993) [28], formula (1) shows that WYt depends on WXt, so the instrumental variable WΔXt can be selected by the first-order difference method for WΔYt. Since ΔXt is a strictly exogenous variable, it is not related to the compound error term ΔEt, so corresponding to the difference model and the level model, the instrumental variables satisfy the following moment conditions:
(11)EWlΔXt′ΔΕt=0,t=3,⋯,T,EWlΔXt′Εt=0,t=3,⋯,T.(3) GMM Estimation. When we estimate the parameters of the spatio-temporal stochastic frontier model, we use the system GMM method similar to the general dynamic panel model to construct the spatial system GMM estimation. Unlike the system GMM, the IVs of the spatial system GMM are composed of time lag variable and spatial lag variable.For each period oft, the moment condition of J≥K+2 can be given. The moment conditions corresponding to the difference model and the level model can be abbreviated as
(12)EHN,ABt′ΔΕt=0,EHN,Lt′Εt=0.The matricesHN,ABt and HN,Lt are expressed as follows
(13)Yt−s=y100000⋯⋯00⋯00y2y1000⋯⋯00⋯0000y3y2y1⋯⋯00⋯0⋮⋮⋮⋮⋮⋮⋯⋯00⋯0000000⋯⋯yt−2yt−3⋯y1,Xt=x3x4⋮xt,ΔYt−s=Δy100000⋯⋯00⋯00Δy2Δy1000⋯⋯00⋯0000Δy3Δy2Δy1⋯⋯00⋯0⋮⋮⋮⋮⋮⋮⋯⋯00⋯0000000⋯⋯Δyt−2Δyt−3⋯Δy1,ΔXt=Δx3Δx4⋮Δxt.Then,HN,ABt=Yt−s,WYt−s,W2Yt−s,⋯,WLYt−s,WΔXt,ΔXt and
(14)HN,Lt=ΔYt−s,WΔYt−s,W2ΔYt−s,⋯,WLΔYt−s,WXt,Xtt=3,⋯,T.That is,HN,ABt and HN,Lt are matrices with instrumental variables as column vectors, and the subscript N means that the matrix depends on the unit number of individuals. Let HN,AB and HN,L be block diagonal matrices composed of block HN,ABt and HN,Lt, respectively. In order to define the GMM (Spatial Blundell Bond, SBB) estimator of the spatial dynamic panel stochastic frontier model, the difference variables and level variables are combined to define the matrix as follows:
(15)Y∗N=ΔY′N,Y′N,Z∗N=ΔZ′N,Z′N,Ε∗N=ΔΕ′N,Ε′N.The instrument variable matrix is(16)HN,BB=diagHN,AB,HN,L=HN,AB00HN,L,where HN,AB is the instrumental matrix of spatial difference GMM estimation, and HN,L is the instrumental matrix of spatial level GMM estimation. The weight matrix is
(17)GN,BB=diagGN,AB,IT−2⊗IN=GN,AB00IT−2⊗IN.This diagonal of the matrix is composed of the weight matrix defined in the process of spatial difference GMM estimation and an identity matrix, whereGN,AB=IN⊗G is NT−2×NT−2 weight matrix which elements are
(18)Gij=2,i=j,−1,i=j+1,−1,j=i+1,0,others.This weight matrix is proposed by Arellano and Bond (1991) [25] which is further define the weight matrix:
(19)AN,BB=H′N,BBGN,BBHN,BBN−1.Through the above process, combining the spatial difference equation with the spatial level equation, we get the spatial system GMM estimation process. Get the objective function of generalized moment estimation for spatial system as follows:(20)1NH′N,BBΔΕ∗N′AN,BBH′N,BBΔΕ∗N=1NH′N,BBΔY∗N−ΔZ∗NθN′AN,BBH′N,BBΔY∗N−ΔZ∗NθN.The one-stage SBB estimator ofθ can be obtained by minimizing Equation (20):
(21)θ^SBB=Z∗′NH′N,BBA∗NH′N,BBZ∗NZ∗′NH′N,BBA∗NH′N,BBY∗N.Equation (21) can also be called the spatial system GMM estimator.(4) Improvement of Instrumental Variable Matrix. The instrumental variable matrix constructed in accordance with the above method has a high dimension and grows exponentially as the values of T and L increase. In order to reduce the dimension of the instrumental variable matrix and avoid overfitting the instrumental variable, we can simplify it by using the “condensing instrumental variable matrix” proposed by Beck and Levine (2004) [29].We still sets=2 and L=1 in the GMM instrumental variable matrix of the space system, and the corresponding condensed instrumental variable matrix is
(22)H′N,BB=H1AB⋯HNABH1L⋯HNL,where HiL (i=1,⋯,N) is the instrumental variable quantum matrix of the level model corresponding to the i individual.
(23)HiL=Δy100⋯WΔy100⋯WΔX3ΔX3Δy2Δy10⋯WΔy2WΔy10⋯WΔX4ΔX4Δy3Δy2Δy1⋯WΔy3WΔy2WΔy1⋯WΔX5ΔX5⋮⋮⋮⋱⋮⋮⋮⋱⋮⋮.


## 2.3. Estimation of Error Term Parameters

### 2.3.1. Estimation ofρ,σv2

After the estimator of the structural parameterθN=λ1,λ2,γ,Β′ in model (1) is obtained in the first stage, the model residual Ε^t=Yt−Ztθ^N can be further obtained, in which the Zt=WYt,WYt−1,Yt−1,Xt is the vector set composed of all explanatory variables in model (1). Consistent GMM estimation can be obtained by using residual Ε^t and modifying the moment condition proposed by Kapoor et al. (2007) [18]. The specific process is as follows.According to the assumptions of the model, the individual effect of the model is the inefficiency term. According to the covariance structure of the compound error term, it can be known that(24)Covεε′=σu2IN⊗JT+σv2INT.Introduce transformation matrix:(25)Q0=IT−JTT⊗IN,Q1=JTT⊗IN,where IT and IN are the identity matrix of order T and order N, JT=eTeT′ is the matrix of order T×T, and the elements of that are all 1. Properties of transformation matrix Q0: (i) trQ0=NT−1; (ii) Q0eT⊗IN=0; (iii) Q0IT⊗RN=IT⊗RNQ0 (where RN is any N×N matrix). From the properties (ii), we can further deduce the special properties (iv) Q0ε=Q0v of matrix Q0 in this paper, where ε and ν are the corresponding error terms in model (1).Let(26)Ε¯=IT⊗MΕ;Ε¯¯=IT⊗MΕ¯;ε¯=IT⊗Mε.Then,(27)ε=Ε−ρΕ¯,ε¯=Ε¯−ρΕ¯¯.Based on the above transformation and referring to the first three of the six moment conditions given by Kapoor et al. (2007) [18] and related properties, the following three moment conditions are given in this paper:
(28)Eε′Q0ε=σv2NT−1,Eε¯′Q0ε¯=σv2NT−1trM′M,Eε¯′Q0ε=0.To further integrate the above moment conditions, we can get that(29)E1NT−1ε′Q0ε1NT−1ε¯′Q0ε¯1NT−1ε¯′Q0ε=σv21Nσv2trM′M0.Substitute Equations (26) and (27) into Equation (29) to obtain that
(30)1NE2T−1Ε¯′Q0E−1T−1Ε¯′Q0Ε¯12T−1Ε¯′Q0Ε¯¯−1T−1Ε¯¯′Q0Ε¯¯trM′M1T−1Ε¯′Q0Ε¯¯+Ε¯′Q0Ε¯−1T−1Ε¯′Q0Ε¯¯0ρρ2σv2=1N1T−1E′Q0E1T−1Ε¯′Q0Ε¯1T−1E′Q0Ε¯.The residualΕ^t=Yt−Ztθ^N estimated in the first stage is substituted into Ε¯ and Ε¯¯ in Equation (30) to obtain the sample moment equation. In the sample moment equation, the estimated value ρ^,σ^v2 of ρ and σv2 can be solved by the following objective function:
(31)ρ^,σ^v2=argminGNρ,ρ2,σv2−gN′GNρ,ρ2,σv2−gN,ρ∈−a,a,σv2∈0,bv,where
(32)GN=1N2T−1Ε∧¯′Q0Ε^−1T−1Ε∧¯′Q0Ε^¯12T−1Ε∧¯′Q0Ε^¯¯−1T−1Ε∧¯¯′Q0Ε^¯¯trM′M1T−1Ε∧¯′Q0Ε^¯¯+Ε∧¯′Q0Ε^¯−1T−1Ε∧¯′Q0Ε^¯¯0,gN=1N1T−1Ε∧′Q0Ε^1T−1Ε∧¯′Q0Ε^¯1T−1Ε∧′Q0Ε^¯.


### 2.3.2. Estimation ofσu2

The fourth moment condition given by Kapoor et al. (2007) [18] is
(33)Eε′Q1ε=NTσu2+Nσv2=Nσ12,where σ12=Tσu2+σv2. However, considering the characteristics of the stochastic frontier model, the compound error term ε obeys the asymmetric distribution of the expected nonzero, so the moment condition (33) cannot be directly applied, and the following formula can be proved:
(34)COVε′Q1ε=NTσu2+Nσv2=Nσ12.The estimator of parameterσ12 can be obtained by combining Equation (34) with Equation (31):
(35)σ^12=1NCOVΕ∧−ρ∧Ε∧¯′Q1Ε^−ρ^Ε^¯=1T−1COVΕ∧′Q1Ε^−2T−1COVΕ∧¯′Q1Ε^·ρ^−1T−1COVΕ∧¯′Q1Ε^¯ρ∧2.Substituteσ^v2 and σ^12 into (33) to obtain the moment estimator σ^u2 of σu2.


## 2.3.1. Estimation ofρ,σv2

After the estimator of the structural parameterθN=λ1,λ2,γ,Β′ in model (1) is obtained in the first stage, the model residual Ε^t=Yt−Ztθ^N can be further obtained, in which the Zt=WYt,WYt−1,Yt−1,Xt is the vector set composed of all explanatory variables in model (1). Consistent GMM estimation can be obtained by using residual Ε^t and modifying the moment condition proposed by Kapoor et al. (2007) [18]. The specific process is as follows.According to the assumptions of the model, the individual effect of the model is the inefficiency term. According to the covariance structure of the compound error term, it can be known that(24)Covεε′=σu2IN⊗JT+σv2INT.Introduce transformation matrix:(25)Q0=IT−JTT⊗IN,Q1=JTT⊗IN,where IT and IN are the identity matrix of order T and order N, JT=eTeT′ is the matrix of order T×T, and the elements of that are all 1. Properties of transformation matrix Q0: (i) trQ0=NT−1; (ii) Q0eT⊗IN=0; (iii) Q0IT⊗RN=IT⊗RNQ0 (where RN is any N×N matrix). From the properties (ii), we can further deduce the special properties (iv) Q0ε=Q0v of matrix Q0 in this paper, where ε and ν are the corresponding error terms in model (1).Let(26)Ε¯=IT⊗MΕ;Ε¯¯=IT⊗MΕ¯;ε¯=IT⊗Mε.Then,(27)ε=Ε−ρΕ¯,ε¯=Ε¯−ρΕ¯¯.Based on the above transformation and referring to the first three of the six moment conditions given by Kapoor et al. (2007) [18] and related properties, the following three moment conditions are given in this paper:
(28)Eε′Q0ε=σv2NT−1,Eε¯′Q0ε¯=σv2NT−1trM′M,Eε¯′Q0ε=0.To further integrate the above moment conditions, we can get that(29)E1NT−1ε′Q0ε1NT−1ε¯′Q0ε¯1NT−1ε¯′Q0ε=σv21Nσv2trM′M0.Substitute Equations (26) and (27) into Equation (29) to obtain that
(30)1NE2T−1Ε¯′Q0E−1T−1Ε¯′Q0Ε¯12T−1Ε¯′Q0Ε¯¯−1T−1Ε¯¯′Q0Ε¯¯trM′M1T−1Ε¯′Q0Ε¯¯+Ε¯′Q0Ε¯−1T−1Ε¯′Q0Ε¯¯0ρρ2σv2=1N1T−1E′Q0E1T−1Ε¯′Q0Ε¯1T−1E′Q0Ε¯.The residualΕ^t=Yt−Ztθ^N estimated in the first stage is substituted into Ε¯ and Ε¯¯ in Equation (30) to obtain the sample moment equation. In the sample moment equation, the estimated value ρ^,σ^v2 of ρ and σv2 can be solved by the following objective function:
(31)ρ^,σ^v2=argminGNρ,ρ2,σv2−gN′GNρ,ρ2,σv2−gN,ρ∈−a,a,σv2∈0,bv,where
(32)GN=1N2T−1Ε∧¯′Q0Ε^−1T−1Ε∧¯′Q0Ε^¯12T−1Ε∧¯′Q0Ε^¯¯−1T−1Ε∧¯¯′Q0Ε^¯¯trM′M1T−1Ε∧¯′Q0Ε^¯¯+Ε∧¯′Q0Ε^¯−1T−1Ε∧¯′Q0Ε^¯¯0,gN=1N1T−1Ε∧′Q0Ε^1T−1Ε∧¯′Q0Ε^¯1T−1Ε∧′Q0Ε^¯.


## 2.3.2. Estimation ofσu2

The fourth moment condition given by Kapoor et al. (2007) [18] is
(33)Eε′Q1ε=NTσu2+Nσv2=Nσ12,where σ12=Tσu2+σv2. However, considering the characteristics of the stochastic frontier model, the compound error term ε obeys the asymmetric distribution of the expected nonzero, so the moment condition (33) cannot be directly applied, and the following formula can be proved:
(34)COVε′Q1ε=NTσu2+Nσv2=Nσ12.The estimator of parameterσ12 can be obtained by combining Equation (34) with Equation (31):
(35)σ^12=1NCOVΕ∧−ρ∧Ε∧¯′Q1Ε^−ρ^Ε^¯=1T−1COVΕ∧′Q1Ε^−2T−1COVΕ∧¯′Q1Ε^·ρ^−1T−1COVΕ∧¯′Q1Ε^¯ρ∧2.Substituteσ^v2 and σ^12 into (33) to obtain the moment estimator σ^u2 of σu2.


## 2.4. Spatial Correction of Estimators

Although it can be proved that the estimator (21) is a consistent estimator, it can also be proved that the consistency of the GMM estimator can be guaranteed even if the model has spatial error autocorrelation. However, the estimator (21) cannot solve the spatial dependence of the error term, and the variance of the estimator is relatively large. After obtaining the consistent estimator of ρ by Equation (31), the consistent estimator can be obtained by a correcting transformation. According to the spatial correction method given by Jacobs et al. (2009) [19], the estimator obtained in the first step was corrected.The estimatorρ^ obtained from Equation (31) is used to construct matrix I−ρ^M, and left the difference GMM and the explained variables and the instrumental variables matrix of the system GMM estimation, if
(36)Y~=I−ρ^MY,ΔY~N=I−ρ^MΔYN,Y~−1=I−ρ^MY−1,W~Y=I−ρ^MWY,X~=I−ρ^MX.The corresponding explanatory variable setZt=WYt,WYt−1,Yt−1,Xt is corrected as
(37)Z~=I−ρ^MZ∗N,ΔZ~=I−ρ^MΔZN.The instrumental variable matrix and weight matrix corresponding to the GMM estimation of the spatial system are corrected as follows:(38)H~SBB=I−ρ^MHN,BB,A~SBB=H~′SBBH~SBB−1.Then, the corrected system GMM estimator is(39)θ~SBB=Z~′SBBH~SBBA~SBBH~′SBBZ~SBBZ~′SBBH~SBBA~SBBH~′SBBY~SBB.


## 3. Results and Discussion

### 3.1. Properties of the Estimators

#### 3.1.1. Properties of the Estimatorθ~SBB

According to the Extremum Consistency Theorem [30] (see Appendix B), the estimators θ^SBB and θ~SBB, obtained by Equations (21) and (27), are consistent.Proof.
see AppendixC.


#### 3.1.2. Properties of the Estimatorsρ^, σ^v2, σ^u2, and θ~SBB

It can be proved that the estimatorsρ^, σ^v2, σ^u2, and θ~SBB are consistent. The proof of the consistency of ρ^, σ^v2, and σ^u2 is similar to that in Kapoor et al. (2007) [18], and it is omitted here. The consistency of θ~SBB can be derived from the consistency of ρ^, σ^v2, and θ^SBB. Similar to the consistency of GLS estimates, the correcting transformation does not affect the consistency of the estimator θ~SBB.


### 3.2. Crossvalidation Scheme and Selection of Spatial Weight Matrix

In order to avoid affecting the accuracy of model estimation due to the choice of spatial weight matrix, the optimal spatial weight matrix was selected by crossvalidation. This is a widely used model selection and generalization method in machine learning. However, since the data used is panel data and the model used is spatio-temporal model, the structural features of spatio-temporal data may be destroyed if the training set and validation set are generated by hold-out method and LOOCV orK-folded crossvalidation. Therefore, this paper considers a stratified crossvalidation approach. For the spatio-temporal data, if N and T are assumed to be the number of spatial individuals and the number of periods contained in the observed samples, respectively, and the rest of the conventions on independent variables and dependent variables are the same as Equation (1), stratified crossvalidation can include the following three forms.


#### 3.2.1. Leave-One-Out Crossvalidation for the Time Dimension (TLOOCV)

Select the datet as the validation set and the rest T−1 of the date as the training set. Let C1,C2,⋯,CT denote, respectively, the index values of the observations contained in period t (t=1,2,⋯,T), and N1,N2,⋯,NT the number of observations contained in period t (t=1,2,⋯,T). Let n1,n2,⋯,nT denote the number of the observations in part t. Do the above for each period t=1,2,⋯,T in turn and calculate
(40)CVT=1T∑t=1TMSEt,where MSEt=∑i=1Ntyit−y∧it2/Nt, and y^it is the fitting value of the ith observed value in period Tyit.This crossvalidation method is suitable when the number of periodsT is not too large.


#### 3.2.2.K-Fold Pooled Crossvalidation for the Spatial Dimension (SK-Fold PCV)

When the total number of periodsT is large and the number of spatial individuals N is also large, this method is suitable for use. All observed values in each period t (t=1,2,⋯,T) were randomly divided into K groups of equal size (i.e., the subsample size of each group was N/K), and a group was randomly selected from each period t (t=1,2,⋯,T) to obtain NT/K observed values combined as the validation set and NK−1T/K remaining observed values in each period combined as the training set. Do the above for each period sequence and calculate
(41)CVK=1K∑k=1KMSEk,where MSEk=∑i=1,t=1Nt,TKyit−y∧it2/TNt, and y^it is the fitting value of the ith observed value in period Tyit.


#### 3.2.3. Leave-One-Out Crossvalidation for the Spatial Dimension (SLOOCV)

When the crossvalidation method presented in Section3.2.2 and the condition K=N attached, it can be called K-fold pooled crossvalidation for the spatial dimension (SN-fold PCV) or stratified leave-one crossvalidation (SLOOCV). When the total number of periods T is large and the number of spatial individuals N is small, this method is suitable for use.


#### 3.2.4. Determination of Weight Matrix and Industrial Efficiency Measure

To discuss industrial efficiency from the perspective of spatial statistics or spatial data mining, a good spatial weight matrix should be determined first. In this paper, the spatial lag production function is chosen as the basic model, and the spatial weight matrix involved in the construction of the model can take various alternative forms. In a data-driven way, the training samples were imported into the model for parameter estimation, and then, the most appropriate spatial weight matrix was determined by stratified crossvalidation. To determine whether the spatial model is selected for analysis, the spatial correlation test is further carried out. If there is a strong spatial correlation, a spatial stochastic frontier model or a spatio-temporal stochastic frontier model will be established; if the spatial correlation is weak, an ordinary panel stochastic frontier model will be selected. After the estimation is completed, the best performing model is used to measure the technical efficiency. The flow chart of the entire analysis process is shown in Figure1.Figure 1
Spatio-temporal industrial efficiency measurement process.


## 3.1. Properties of the Estimators

### 3.1.1. Properties of the Estimatorθ~SBB

According to the Extremum Consistency Theorem [30] (see Appendix B), the estimators θ^SBB and θ~SBB, obtained by Equations (21) and (27), are consistent.Proof.
see AppendixC.


### 3.1.2. Properties of the Estimatorsρ^, σ^v2, σ^u2, and θ~SBB

It can be proved that the estimatorsρ^, σ^v2, σ^u2, and θ~SBB are consistent. The proof of the consistency of ρ^, σ^v2, and σ^u2 is similar to that in Kapoor et al. (2007) [18], and it is omitted here. The consistency of θ~SBB can be derived from the consistency of ρ^, σ^v2, and θ^SBB. Similar to the consistency of GLS estimates, the correcting transformation does not affect the consistency of the estimator θ~SBB.


## 3.1.1. Properties of the Estimatorθ~SBB

According to the Extremum Consistency Theorem [30] (see Appendix B), the estimators θ^SBB and θ~SBB, obtained by Equations (21) and (27), are consistent.Proof.
see AppendixC.


## 3.1.2. Properties of the Estimatorsρ^, σ^v2, σ^u2, and θ~SBB

It can be proved that the estimatorsρ^, σ^v2, σ^u2, and θ~SBB are consistent. The proof of the consistency of ρ^, σ^v2, and σ^u2 is similar to that in Kapoor et al. (2007) [18], and it is omitted here. The consistency of θ~SBB can be derived from the consistency of ρ^, σ^v2, and θ^SBB. Similar to the consistency of GLS estimates, the correcting transformation does not affect the consistency of the estimator θ~SBB.


## 3.2. Crossvalidation Scheme and Selection of Spatial Weight Matrix

In order to avoid affecting the accuracy of model estimation due to the choice of spatial weight matrix, the optimal spatial weight matrix was selected by crossvalidation. This is a widely used model selection and generalization method in machine learning. However, since the data used is panel data and the model used is spatio-temporal model, the structural features of spatio-temporal data may be destroyed if the training set and validation set are generated by hold-out method and LOOCV orK-folded crossvalidation. Therefore, this paper considers a stratified crossvalidation approach. For the spatio-temporal data, if N and T are assumed to be the number of spatial individuals and the number of periods contained in the observed samples, respectively, and the rest of the conventions on independent variables and dependent variables are the same as Equation (1), stratified crossvalidation can include the following three forms.


### 3.2.1. Leave-One-Out Crossvalidation for the Time Dimension (TLOOCV)

Select the datet as the validation set and the rest T−1 of the date as the training set. Let C1,C2,⋯,CT denote, respectively, the index values of the observations contained in period t (t=1,2,⋯,T), and N1,N2,⋯,NT the number of observations contained in period t (t=1,2,⋯,T). Let n1,n2,⋯,nT denote the number of the observations in part t. Do the above for each period t=1,2,⋯,T in turn and calculate
(40)CVT=1T∑t=1TMSEt,where MSEt=∑i=1Ntyit−y∧it2/Nt, and y^it is the fitting value of the ith observed value in period Tyit.This crossvalidation method is suitable when the number of periodsT is not too large.


### 3.2.2.K-Fold Pooled Crossvalidation for the Spatial Dimension (SK-Fold PCV)

When the total number of periodsT is large and the number of spatial individuals N is also large, this method is suitable for use. All observed values in each period t (t=1,2,⋯,T) were randomly divided into K groups of equal size (i.e., the subsample size of each group was N/K), and a group was randomly selected from each period t (t=1,2,⋯,T) to obtain NT/K observed values combined as the validation set and NK−1T/K remaining observed values in each period combined as the training set. Do the above for each period sequence and calculate
(41)CVK=1K∑k=1KMSEk,where MSEk=∑i=1,t=1Nt,TKyit−y∧it2/TNt, and y^it is the fitting value of the ith observed value in period Tyit.


### 3.2.3. Leave-One-Out Crossvalidation for the Spatial Dimension (SLOOCV)

When the crossvalidation method presented in Section3.2.2 and the condition K=N attached, it can be called K-fold pooled crossvalidation for the spatial dimension (SN-fold PCV) or stratified leave-one crossvalidation (SLOOCV). When the total number of periods T is large and the number of spatial individuals N is small, this method is suitable for use.


### 3.2.4. Determination of Weight Matrix and Industrial Efficiency Measure

To discuss industrial efficiency from the perspective of spatial statistics or spatial data mining, a good spatial weight matrix should be determined first. In this paper, the spatial lag production function is chosen as the basic model, and the spatial weight matrix involved in the construction of the model can take various alternative forms. In a data-driven way, the training samples were imported into the model for parameter estimation, and then, the most appropriate spatial weight matrix was determined by stratified crossvalidation. To determine whether the spatial model is selected for analysis, the spatial correlation test is further carried out. If there is a strong spatial correlation, a spatial stochastic frontier model or a spatio-temporal stochastic frontier model will be established; if the spatial correlation is weak, an ordinary panel stochastic frontier model will be selected. After the estimation is completed, the best performing model is used to measure the technical efficiency. The flow chart of the entire analysis process is shown in Figure1.Figure 1
Spatio-temporal industrial efficiency measurement process.


## 3.2.1. Leave-One-Out Crossvalidation for the Time Dimension (TLOOCV)

Select the datet as the validation set and the rest T−1 of the date as the training set. Let C1,C2,⋯,CT denote, respectively, the index values of the observations contained in period t (t=1,2,⋯,T), and N1,N2,⋯,NT the number of observations contained in period t (t=1,2,⋯,T). Let n1,n2,⋯,nT denote the number of the observations in part t. Do the above for each period t=1,2,⋯,T in turn and calculate
(40)CVT=1T∑t=1TMSEt,where MSEt=∑i=1Ntyit−y∧it2/Nt, and y^it is the fitting value of the ith observed value in period Tyit.This crossvalidation method is suitable when the number of periodsT is not too large.


## 3.2.2.K-Fold Pooled Crossvalidation for the Spatial Dimension (SK-Fold PCV)

When the total number of periodsT is large and the number of spatial individuals N is also large, this method is suitable for use. All observed values in each period t (t=1,2,⋯,T) were randomly divided into K groups of equal size (i.e., the subsample size of each group was N/K), and a group was randomly selected from each period t (t=1,2,⋯,T) to obtain NT/K observed values combined as the validation set and NK−1T/K remaining observed values in each period combined as the training set. Do the above for each period sequence and calculate
(41)CVK=1K∑k=1KMSEk,where MSEk=∑i=1,t=1Nt,TKyit−y∧it2/TNt, and y^it is the fitting value of the ith observed value in period Tyit.


## 3.2.3. Leave-One-Out Crossvalidation for the Spatial Dimension (SLOOCV)

When the crossvalidation method presented in Section3.2.2 and the condition K=N attached, it can be called K-fold pooled crossvalidation for the spatial dimension (SN-fold PCV) or stratified leave-one crossvalidation (SLOOCV). When the total number of periods T is large and the number of spatial individuals N is small, this method is suitable for use.


## 3.2.4. Determination of Weight Matrix and Industrial Efficiency Measure

To discuss industrial efficiency from the perspective of spatial statistics or spatial data mining, a good spatial weight matrix should be determined first. In this paper, the spatial lag production function is chosen as the basic model, and the spatial weight matrix involved in the construction of the model can take various alternative forms. In a data-driven way, the training samples were imported into the model for parameter estimation, and then, the most appropriate spatial weight matrix was determined by stratified crossvalidation. To determine whether the spatial model is selected for analysis, the spatial correlation test is further carried out. If there is a strong spatial correlation, a spatial stochastic frontier model or a spatio-temporal stochastic frontier model will be established; if the spatial correlation is weak, an ordinary panel stochastic frontier model will be selected. After the estimation is completed, the best performing model is used to measure the technical efficiency. The flow chart of the entire analysis process is shown in Figure1.Figure 1
Spatio-temporal industrial efficiency measurement process.


## 4. The Efficiency of the High and New Technology Industries in China

China has developed high and new technology industries for many years in order to transform the economic growth mode, cultivating knowledge and technology intensive new companies with great growth potential and low resources consumption that provide a sustainable development. As the technology of such industries disseminates, some issues emerge, such as spatial technology spillover, continuity of technological upgrade, and delay from research and development to market acceptance. In this paper, we analyze the efficiency of this kind of industries by the above analysis process.


### 4.1. Introduction to the Model and Data

In the framework of spatio-temporal model analysis, the spatial lag production function model based on the Cobb-Douglas production function is chosen as the basic model for weight selection. The matrix form of the model is as follows:(42)lnY=ρWNlnY+lnk,lnlΒ+u,where lnY, lnk, and lnl are, respectively, the logarithm vectors of the main business income, the assets investment, and the mean number of employees of the high and new technology industries in every province in China; WN is the spatial weight matrix of the 31 provinces in China. Among many spatial weight matrices, we choose the three most commonly used spatial weight matrices in economic problems and their combination forms. Various weight matrices and their interpretations are shown in Table 1.Table 1
Various weight matrices to be selected.
ExpressionMatrix nameMeaningW1Position adjacency weight matrixThe weight matrix is constructed by rook, bishop, and queen position adjacence, and queen adjacence matrix is selected in this paper.Matrix element wij=0Regioniandjarenotadjacent1RegioniandjareadjacentW2Geographical distance weight matrixThe weight matrix is constructed by geographical distance between regions, and this paper constructs the weight matrix by reciprocal distance between the centers of provincial capitals in China.Matrix element wij=1/dij, where dij is the geographical distance between regions i and jW3Economic distance weight matrixThe weight matrix is constructed by the difference of economic level among different regions, that is, the smaller the economic gap, the stronger the spatial correlation. In this paper, the GDP of each province in China is used as the proxy variable of economic development level to construct the matrixMatrix element wij=1/GDP¯i−GDP¯ji≠j0i=j (GDP¯i represents the annual average GDP of region i)W4Adjacency and distance combination weight matrixMatrix elementwij=0Regioniandjarenotadjacent1/dijRegioniandjareadjacentW5Adjacency and economic combination weight matrixMatrix elementwij=0Regioniandjarenotadjacent1/GDP¯i−GDP¯jRegioniandjareadjacentW6Weight matrix of distance and economy combinationMatrix elementwij=1/dij·GDP¯i/GDP¯j (where dij is the geographical distance between region i and j; GDP¯i represents the annual average GDP of region i)Before the fitting, each weight matrix was row-standardized, and the optimal weight matrix determined by crossvalidation was implemented to establish the stochastic frontier model. Starting from the ordinary panel stochastic frontier model and considering the spatial correlations, we construct the static panel spatial stochastic frontier model and the spatiotemporal stochastic frontier model. In order to determine if the model variables present spatial correlation, we let them go through a spatial correlation test, and to determine if there should be a time lag term in the model, we test the significance. Comparing the results of the three models and selecting the one that provides the best fit, we estimate the technology efficiency of the high and new technology industries in every province.The matrix forms of the three models are as follows:Ordinary panel stochastic frontier model:(43)lnY=lnk,lnlΒ+v−u.Static panel space stochastic frontier model:(44)lnY=λIT⊗WNlnY+lnk,lnlΒ+v−u,v=ρIT⊗MNv+ξ.Spatiotemporal stochastic frontier model:(45)lnY=λIT⊗WlnY+πIT⊗WlnY‐1+lnk,lnlΒ+Ε,v=ρIT⊗MNv+ξ,where v is a general vector of stochastic error, u is the inefficiency term, IT is the identity matrix of order T, λ and ρ are the spatial correlation coefficients of the corresponding equations, π is the spatiotemporal time lag coefficient of the spatiotemporal stochastic frontier model, and B is the regression coefficient vector.The development plan of high and new technology industry in China started in 1988, but due to the relatively slow progress in the beginning, the scale development of this industry did not start until the beginning of the twenty-first century. For this reason, we have chosen as research sample the panel data of the high and new technology industries of the 31 provinces in China from 2001 to 2018. The data of capital and labor input factors have been taken from “China high-tech industry yearbook.” Descriptive statistics are shown in Table2.Table 2
Descriptive statistics of output value and factor input index of high-tech industry.
VariableAverageStandard deviationMedianMaximumY30.28353.71213.615237.170k31.68254.75715.884241.980l39.47775.53018.817372.000Figure2 is the histogram drawn by taking the intragroup mean of data in each region according to year.Figure 2
Histogram of the mean value of output value and factor input of high-tech industries in each province.Figure2 shows the difference in investment and average development level of high-tech industries in different provinces of China from 2001 to 2018. As can be seen from the figure, Guangdong, Jiangsu, and Shandong provinces have the highest input and output levels of high-tech industries, and the differences among these three provinces are also very large. In the past 18 years, the average output value of the high-tech industry of Guangdong, which ranks first, reached 23.7 billion yuan, while that of Shandong, which ranks third, reached only 6.997 billion yuan, less than one third of that of Guangdong. In terms of the development and distribution of high-tech industries nationwide, the gap between provinces is even more obvious. The average output value of the high-tech industries in Tibet, which ranks the last, is only 0.085 billion yuan, less than 1/1000 of that of Guangdong.


### 4.2. Crossvalidation Results and the Selection of Weight Matrix

According to the characteristics of the data obtained, the time limit contained in the data is 18 years, which is relatively short and smaller than the number of regions. Therefore, the leave-one-out crossvalidation for the time dimension (TLOOCV) method was chosen. Each weight matrix in Table1 was introduced into model (31), the training set data were imported into the model one by one for fitting, and then, Equation (40) was calculated to obtain the CV statistics corresponding to each weight matrix. The calculation results are shown in Table 3.Table 3
CV statistics corresponding to various weight matrices.
Weight matricesW1W2W3W4W5W6CV statistics526.932403.9296333.7308157.5527320.589104.8576By comparing the calculation results of CV statistics of validation set, it can be found to be the optimal spatial weight matrix required by this paper.


### 4.3. Empirical Results

To ensure that spatial econometrics is applicable to the problem we are studying, we need to test the spatial correlation of the variables we are interested in. The most popular method to measure spatial autocorrelation is Moran’s indexI (Moran’s I):I=∑i=1n∑j=1nwijxi−x¯xj−x¯/S2∑i=1n∑j=1nwij, where S2 is the sample variance, wij is the (i,j) element of the spatial weight matrix (used to measure the distance between region i and region j), and ∑i=1n∑j=1nwij is the sum of all spatial weights.The value of Moran’sI is generally between -1 and 1, and its greater than 0 indicates positive autocorrelation. That is, the high value is adjacent to the high value and the low value is adjacent to the low value. Less than 0 means negative autocorrelation. That is, a high value is adjacent to a low value. If the Moran’s I is close to 0, then the spatial distribution is random, and there is no spatial autocorrelation.To test the existence of spatial correlations in the variables of the high and new technology industry, we calculate the global Moran’sI indices of the production values of the industry from 2011 to 2018 (Table 4).Table 4
Moran’sI indices among the production values of the high and new technology industry.
IndexP valueIndexP valueIndexP valueIndexP valueYear2011201220132014MoranI0.286P≤0.010.321P≤0.0010.340P≤0.0010.302P≤0.01Year2015201620172018MoranI0.312P≤0.010.298P≤0.010.305P≤0.010.317P≤0.01From the results of Moran’sI index calculation, we found that the P value of the index is smaller than 0.01 for every year, demonstrating that the index is significant below 1% for every year, and the average Moran’s I index is also significant for every year. The Moran’s I index reached the minimum value 0.286 in 2011 and the maximum value 0.340 in 2013. We observe that the production value of the high and new technology industries of every province shows significant spatial correlation for every year and conclude that the production values of these industries of the provinces in China have apparent spatial aggregating effects. Furthermore, the testing of spatial correlation and the location quotient calculation both demonstrated that the high and new technology industries in different regions of China have apparent spatial correlation. We therefore choose the panel space stochastic frontier model for the analysis.Table5 presents the estimation results of the static, spatiotemporal stochastic frontier model, where the static model was further analyzed by considering fixed and stochastic effects.Table 5
Estimation results by the static panel space and spatiotemporal stochastic frontier model.
Dependent variableLnyModelStatic panel spatial stochastic frontierSpatiotemporal stochastic frontierVariableFixed effectsRandom effectsCoefficientsP valueCoefficientsP valueCoefficientsP valuelnk0.397∗∗∗0.0010.401∗∗0.0010.405∗∗0.008lnl0.686∗∗∗0.0080.713∗∗0.0170.718∗∗∗0.005λ0.435∗∗0.0120.576∗∗∗0.0020.369∗∗∗0.005π————0.205∗∗0.011ρ0.551∗∗∗0.0070.698∗∗0.0130.585∗∗0.017σu22.65E-050.3370.656σv20.3870.1320.180γ6.850E-040.7190.786Hausmanstatistics-39282.77D statistics10.668∗∗In the table,∗∗ indicates those that can pass the test with 5% significance level.The spatial autoregressive coefficientsλ and ρ of the three models can all pass the significance test. From this and the spatial correlation test, one can conclude that the spatial stochastic frontier model is more reasonable. The estimated spatial autoregressive coefficients of the three models are all positive, which implied that the spatial effects have a positive impact on the development of the high and new technology industries. The negative value of the Hausman statistics of the static panel space model implies that the random effect model should be chosen. The random effect σu2 value of the static panel space stochastic frontier model is far greater than the σv2 value, and the value of γ is 0.719, implying that there apparently exists technical inefficiency. Comparing the spatiotemporal stochastic frontier model with the estimation of the static panel spatial random effect, the spatiotemporal lag coefficient of the former can pass the 10% significance test, obtaining that the spatiotemporal lag term in the model has a significant function. The distance statistic of the spatiotemporal stochastic frontier model pass the 5% significance distance test, proving that the spatiotemporal stochastic frontier model is globally significant. Moreover, the estimation of the γ value of the spatiotemporal stochastic frontier model is higher than that of the static panel spatial stochastic frontier model, demonstrating that the inefficiency term of the spatiotemporal stochastic frontier model has a more significant function. All the impact factor variables of the technical inefficiency terms of the analysis of the two models can at least pass the 5% significance test, and the signs of the regression coefficient of the two models are consistent, the numerical values are relatively close. Taken together, the above results all demonstrate that the analysis of the spatiotemporal stochastic frontier model is more reasonable and the development of the high and new technology industry has positive correlations in space and time. Due to space limitation, this paper does not report the annual technical efficiency of the high-tech industries of each province. The spatial-temporal stochastic frontier model is used to calculate the average technical efficiency of the high-tech industries of each province from 2001 to 2018 as follows.It can be seen from the calculation results in Table6 that the average technical efficiency value of the high-tech industries in all provinces in China is less than 1, which indicates that the actual output of the high-tech industries in all provinces has not reached the most effective output level, and there is technological inefficiency in production. The five-year national average technical efficiency level was 0.837, and there are obvious regional differences in the technical efficiency values presented in Table 4. Nine provinces (Beijing, Chongqing, Fujian, Gansu, Jiangsu, Liaoning, Shanghai, Tianjin, and Zhejiang) achieved an average technical efficiency of more than 0.9, seven provinces are located in the eastern region, one in the central region, and one in the western region. There are 11 provinces with average technical efficiency below 0.8, namely, Guangdong, Guangxi, Guizhou, Hebei, Henan, Jiangxi, Ningxia, Qinghai, Shaanxi, Xinjiang, and Tibet. Only one of the provinces is in the east, six in the central region, and four in the west.Table 6
Average technological efficiency of high-tech industries in 31 provinces.
ProvinceTEProvinceTEProvinceTEProvinceTEAnhui0.817Guizhou0.818Hunan0.812Ningxia0.786Beijing0.985Hainan0.876Jilin0.835Qinghai0.733Chongqing0.912Hebei0.786Jiangsu0.951Sichuan0.851Fujian0.936Heilongjiang0.835Jiangxi0.792Shandong0.876Gansu0.932Henan0.752Liaoning0.931Shanghai0.966Guangdong0.755Hubei0.851Neimenggu0.872Shǎnxi0.795Guangxi0.785Shānxi0.787Tianjin0.975Xinjiang0.651Tibet0.668Yunnan0.815Zhejiang0.952Data source: calculated based on pattern II stochastic effects model of a function.


## 4.1. Introduction to the Model and Data

In the framework of spatio-temporal model analysis, the spatial lag production function model based on the Cobb-Douglas production function is chosen as the basic model for weight selection. The matrix form of the model is as follows:(42)lnY=ρWNlnY+lnk,lnlΒ+u,where lnY, lnk, and lnl are, respectively, the logarithm vectors of the main business income, the assets investment, and the mean number of employees of the high and new technology industries in every province in China; WN is the spatial weight matrix of the 31 provinces in China. Among many spatial weight matrices, we choose the three most commonly used spatial weight matrices in economic problems and their combination forms. Various weight matrices and their interpretations are shown in Table 1.Table 1
Various weight matrices to be selected.
ExpressionMatrix nameMeaningW1Position adjacency weight matrixThe weight matrix is constructed by rook, bishop, and queen position adjacence, and queen adjacence matrix is selected in this paper.Matrix element wij=0Regioniandjarenotadjacent1RegioniandjareadjacentW2Geographical distance weight matrixThe weight matrix is constructed by geographical distance between regions, and this paper constructs the weight matrix by reciprocal distance between the centers of provincial capitals in China.Matrix element wij=1/dij, where dij is the geographical distance between regions i and jW3Economic distance weight matrixThe weight matrix is constructed by the difference of economic level among different regions, that is, the smaller the economic gap, the stronger the spatial correlation. In this paper, the GDP of each province in China is used as the proxy variable of economic development level to construct the matrixMatrix element wij=1/GDP¯i−GDP¯ji≠j0i=j (GDP¯i represents the annual average GDP of region i)W4Adjacency and distance combination weight matrixMatrix elementwij=0Regioniandjarenotadjacent1/dijRegioniandjareadjacentW5Adjacency and economic combination weight matrixMatrix elementwij=0Regioniandjarenotadjacent1/GDP¯i−GDP¯jRegioniandjareadjacentW6Weight matrix of distance and economy combinationMatrix elementwij=1/dij·GDP¯i/GDP¯j (where dij is the geographical distance between region i and j; GDP¯i represents the annual average GDP of region i)Before the fitting, each weight matrix was row-standardized, and the optimal weight matrix determined by crossvalidation was implemented to establish the stochastic frontier model. Starting from the ordinary panel stochastic frontier model and considering the spatial correlations, we construct the static panel spatial stochastic frontier model and the spatiotemporal stochastic frontier model. In order to determine if the model variables present spatial correlation, we let them go through a spatial correlation test, and to determine if there should be a time lag term in the model, we test the significance. Comparing the results of the three models and selecting the one that provides the best fit, we estimate the technology efficiency of the high and new technology industries in every province.The matrix forms of the three models are as follows:Ordinary panel stochastic frontier model:(43)lnY=lnk,lnlΒ+v−u.Static panel space stochastic frontier model:(44)lnY=λIT⊗WNlnY+lnk,lnlΒ+v−u,v=ρIT⊗MNv+ξ.Spatiotemporal stochastic frontier model:(45)lnY=λIT⊗WlnY+πIT⊗WlnY‐1+lnk,lnlΒ+Ε,v=ρIT⊗MNv+ξ,where v is a general vector of stochastic error, u is the inefficiency term, IT is the identity matrix of order T, λ and ρ are the spatial correlation coefficients of the corresponding equations, π is the spatiotemporal time lag coefficient of the spatiotemporal stochastic frontier model, and B is the regression coefficient vector.The development plan of high and new technology industry in China started in 1988, but due to the relatively slow progress in the beginning, the scale development of this industry did not start until the beginning of the twenty-first century. For this reason, we have chosen as research sample the panel data of the high and new technology industries of the 31 provinces in China from 2001 to 2018. The data of capital and labor input factors have been taken from “China high-tech industry yearbook.” Descriptive statistics are shown in Table2.Table 2
Descriptive statistics of output value and factor input index of high-tech industry.
VariableAverageStandard deviationMedianMaximumY30.28353.71213.615237.170k31.68254.75715.884241.980l39.47775.53018.817372.000Figure2 is the histogram drawn by taking the intragroup mean of data in each region according to year.Figure 2
Histogram of the mean value of output value and factor input of high-tech industries in each province.Figure2 shows the difference in investment and average development level of high-tech industries in different provinces of China from 2001 to 2018. As can be seen from the figure, Guangdong, Jiangsu, and Shandong provinces have the highest input and output levels of high-tech industries, and the differences among these three provinces are also very large. In the past 18 years, the average output value of the high-tech industry of Guangdong, which ranks first, reached 23.7 billion yuan, while that of Shandong, which ranks third, reached only 6.997 billion yuan, less than one third of that of Guangdong. In terms of the development and distribution of high-tech industries nationwide, the gap between provinces is even more obvious. The average output value of the high-tech industries in Tibet, which ranks the last, is only 0.085 billion yuan, less than 1/1000 of that of Guangdong.


## 4.2. Crossvalidation Results and the Selection of Weight Matrix

According to the characteristics of the data obtained, the time limit contained in the data is 18 years, which is relatively short and smaller than the number of regions. Therefore, the leave-one-out crossvalidation for the time dimension (TLOOCV) method was chosen. Each weight matrix in Table1 was introduced into model (31), the training set data were imported into the model one by one for fitting, and then, Equation (40) was calculated to obtain the CV statistics corresponding to each weight matrix. The calculation results are shown in Table 3.Table 3
CV statistics corresponding to various weight matrices.
Weight matricesW1W2W3W4W5W6CV statistics526.932403.9296333.7308157.5527320.589104.8576By comparing the calculation results of CV statistics of validation set, it can be found to be the optimal spatial weight matrix required by this paper.


## 4.3. Empirical Results

To ensure that spatial econometrics is applicable to the problem we are studying, we need to test the spatial correlation of the variables we are interested in. The most popular method to measure spatial autocorrelation is Moran’s indexI (Moran’s I):I=∑i=1n∑j=1nwijxi−x¯xj−x¯/S2∑i=1n∑j=1nwij, where S2 is the sample variance, wij is the (i,j) element of the spatial weight matrix (used to measure the distance between region i and region j), and ∑i=1n∑j=1nwij is the sum of all spatial weights.The value of Moran’sI is generally between -1 and 1, and its greater than 0 indicates positive autocorrelation. That is, the high value is adjacent to the high value and the low value is adjacent to the low value. Less than 0 means negative autocorrelation. That is, a high value is adjacent to a low value. If the Moran’s I is close to 0, then the spatial distribution is random, and there is no spatial autocorrelation.To test the existence of spatial correlations in the variables of the high and new technology industry, we calculate the global Moran’sI indices of the production values of the industry from 2011 to 2018 (Table 4).Table 4
Moran’sI indices among the production values of the high and new technology industry.
IndexP valueIndexP valueIndexP valueIndexP valueYear2011201220132014MoranI0.286P≤0.010.321P≤0.0010.340P≤0.0010.302P≤0.01Year2015201620172018MoranI0.312P≤0.010.298P≤0.010.305P≤0.010.317P≤0.01From the results of Moran’sI index calculation, we found that the P value of the index is smaller than 0.01 for every year, demonstrating that the index is significant below 1% for every year, and the average Moran’s I index is also significant for every year. The Moran’s I index reached the minimum value 0.286 in 2011 and the maximum value 0.340 in 2013. We observe that the production value of the high and new technology industries of every province shows significant spatial correlation for every year and conclude that the production values of these industries of the provinces in China have apparent spatial aggregating effects. Furthermore, the testing of spatial correlation and the location quotient calculation both demonstrated that the high and new technology industries in different regions of China have apparent spatial correlation. We therefore choose the panel space stochastic frontier model for the analysis.Table5 presents the estimation results of the static, spatiotemporal stochastic frontier model, where the static model was further analyzed by considering fixed and stochastic effects.Table 5
Estimation results by the static panel space and spatiotemporal stochastic frontier model.
Dependent variableLnyModelStatic panel spatial stochastic frontierSpatiotemporal stochastic frontierVariableFixed effectsRandom effectsCoefficientsP valueCoefficientsP valueCoefficientsP valuelnk0.397∗∗∗0.0010.401∗∗0.0010.405∗∗0.008lnl0.686∗∗∗0.0080.713∗∗0.0170.718∗∗∗0.005λ0.435∗∗0.0120.576∗∗∗0.0020.369∗∗∗0.005π————0.205∗∗0.011ρ0.551∗∗∗0.0070.698∗∗0.0130.585∗∗0.017σu22.65E-050.3370.656σv20.3870.1320.180γ6.850E-040.7190.786Hausmanstatistics-39282.77D statistics10.668∗∗In the table,∗∗ indicates those that can pass the test with 5% significance level.The spatial autoregressive coefficientsλ and ρ of the three models can all pass the significance test. From this and the spatial correlation test, one can conclude that the spatial stochastic frontier model is more reasonable. The estimated spatial autoregressive coefficients of the three models are all positive, which implied that the spatial effects have a positive impact on the development of the high and new technology industries. The negative value of the Hausman statistics of the static panel space model implies that the random effect model should be chosen. The random effect σu2 value of the static panel space stochastic frontier model is far greater than the σv2 value, and the value of γ is 0.719, implying that there apparently exists technical inefficiency. Comparing the spatiotemporal stochastic frontier model with the estimation of the static panel spatial random effect, the spatiotemporal lag coefficient of the former can pass the 10% significance test, obtaining that the spatiotemporal lag term in the model has a significant function. The distance statistic of the spatiotemporal stochastic frontier model pass the 5% significance distance test, proving that the spatiotemporal stochastic frontier model is globally significant. Moreover, the estimation of the γ value of the spatiotemporal stochastic frontier model is higher than that of the static panel spatial stochastic frontier model, demonstrating that the inefficiency term of the spatiotemporal stochastic frontier model has a more significant function. All the impact factor variables of the technical inefficiency terms of the analysis of the two models can at least pass the 5% significance test, and the signs of the regression coefficient of the two models are consistent, the numerical values are relatively close. Taken together, the above results all demonstrate that the analysis of the spatiotemporal stochastic frontier model is more reasonable and the development of the high and new technology industry has positive correlations in space and time. Due to space limitation, this paper does not report the annual technical efficiency of the high-tech industries of each province. The spatial-temporal stochastic frontier model is used to calculate the average technical efficiency of the high-tech industries of each province from 2001 to 2018 as follows.It can be seen from the calculation results in Table6 that the average technical efficiency value of the high-tech industries in all provinces in China is less than 1, which indicates that the actual output of the high-tech industries in all provinces has not reached the most effective output level, and there is technological inefficiency in production. The five-year national average technical efficiency level was 0.837, and there are obvious regional differences in the technical efficiency values presented in Table 4. Nine provinces (Beijing, Chongqing, Fujian, Gansu, Jiangsu, Liaoning, Shanghai, Tianjin, and Zhejiang) achieved an average technical efficiency of more than 0.9, seven provinces are located in the eastern region, one in the central region, and one in the western region. There are 11 provinces with average technical efficiency below 0.8, namely, Guangdong, Guangxi, Guizhou, Hebei, Henan, Jiangxi, Ningxia, Qinghai, Shaanxi, Xinjiang, and Tibet. Only one of the provinces is in the east, six in the central region, and four in the west.Table 6
Average technological efficiency of high-tech industries in 31 provinces.
ProvinceTEProvinceTEProvinceTEProvinceTEAnhui0.817Guizhou0.818Hunan0.812Ningxia0.786Beijing0.985Hainan0.876Jilin0.835Qinghai0.733Chongqing0.912Hebei0.786Jiangsu0.951Sichuan0.851Fujian0.936Heilongjiang0.835Jiangxi0.792Shandong0.876Gansu0.932Henan0.752Liaoning0.931Shanghai0.966Guangdong0.755Hubei0.851Neimenggu0.872Shǎnxi0.795Guangxi0.785Shānxi0.787Tianjin0.975Xinjiang0.651Tibet0.668Yunnan0.815Zhejiang0.952Data source: calculated based on pattern II stochastic effects model of a function.


## 5. Conclusion

In this paper, taking into account that the variables to be explained might be affected by the time lag term and the space-time interaction, we develop a dynamic model within the framework of the panel spatial stochastic frontier model. Due to the apparent endogeneity of the model, we use the systematic GMM method to estimate the parameters, choose suitable tool variables according to the model assumptions and variable characteristics, and construct the suitable spatiotemporal stochastic frontier model. We use the extreme value consistence theorem and the uniform law of large numbers (ULLN) to prove the consistency of the structural parameter estimators and of the estimators of the error term distribution parameters. Aiming at the selection of spatial weight matrix of spatio-temporal model, a stratified crossvalidation method is designed to select the most appropriate spatial weight matrix in a data-driven way according to the characteristics of spatio-temporal data. Although the spatial weight matrix selected by supervised learning may not be suitable for analyzing all problems, this data-driven model selection method is undoubtedly valuable and efficient.From the analysis of the stochastic frontier model of high-tech industries in China and the measurement of their technical efficiency, we can draw the following conclusions.There is a spatial positive correlation in the development of high-tech industries between different regions of China. The positive correlation between the output values of these industries in different regions has been obtained by calculating the Global Moran’sI index in each year. The estimation of the spatial panel stochastic frontier model also indicates that the spatial autoregressive coefficient is positive, proving the existence of such a positive correlation which has a positive impact on the development of high-tech industries. There is also a spatial agglomeration effect and a spatial and temporal lag effect in these industries, illustrating that both static spill over and dynamic continuity occur in the development of the high-tech industries in China. The technical efficiency of high-tech industries is relatively low. The strategic emerging industries started earlier in eastern region, but developed more slowly than in the central and western regions.The Chinese economy is at a critical stage of replacing old drivers of growth with new ones and transforming and upgrading industries. The new round of technological and industrial revolution 5.0 has given rise to new technologies, new industries, new forms of business, and new models. In this study, the data mining algorithm based on stochastic frontier is used to calculate industrial efficiency, which is not only suitable for high-tech industry but also helpful to further enrich the research on the efficiency of new industry and new mode and has certain practical significance to promote the steady development of the new round of scientific and technological revolution of industry 5.0.



---

*Source: 7439744-2021-11-08.xml*