A spatial stochastic frontier model including both frontier and error-based spatial cross-sectional dependence

ABSTRACT Spatial dependence in stochastic frontier models is usually handled by modelling the frontier function or the inefficiency error term through the introduction of some spatial components. The model proposed in this paper (SDF-CSD) combines the two different modelling approaches, obtaining a full comprehensive specification that introduces four different sources of spatial cross-sectional dependence. The most appealing feature of the model is that it allows capturing global and local spatial spillover effects while controlling for spatial correlation related to firms’ efficiency and to unobserved but spatially correlated variables. Moreover, it can be estimated using maximum likelihood techniques. Finally, some Monte Carlo simulations were run to test the final sample properties of the new spatial estimator and an application to the Italian agricultural sector is provided.


INTRODUCTION
In the last two decades, scholars started to expand classical stochastic frontier (SF) models in order to include some spatial components. Indeed, firms tend to concentrate in clusters, taking advantage of positive agglomeration externalities due to cooperation, shared ideas and emulation, resulting in increased productivity levels (Baptista, 2000). Thus, producers cannot be regarded as isolated entities and the hypothesis of cross-sectional independence underlying the basic SF model must no longer be considered valid. Until now scholars have introduced spatial dependence in SF models following two different paths: evaluating global and local spatial spillovers affecting the frontier function or considering spatial cross-sectional correlation in the inefficiency and/or in the error terms.
Considering the first approach, Adetutu et al. (2015) are the first to introduce local spatial dependence in SF models through the spatial lag of the exogenous variables. Afterwards, Glass et al. (2016) proposed a spatial Durbin SF specification that accounts for both local and global spatial dependence, including the spatial lag of the input variables (SLX term) and of the dependent variable (SAR term). While the former allows capturing input-sharing and labour-matching phenomena, the latter involves spatial dependence arising from joint practices, knowledge transmission, adoption of common technologies, creation of a collective identity and market-oriented externalities. Further extensions of the spatial Durbin SF model introduced by Glass et al. (2016) were developed by Ramajo and Hewings (2018) and Gude et al. (2018) considering, respectively, a time-varying decay specification for the inefficiency error term and timevarying exogenous influences on both the degree of global and local spatial dependence for each time period. Finally, Tsukamoto (2019), following the specification of Battese and Coelli (1995) for non-spatial SF models, introduced the possibility to consider some determinants of firms' technical inefficiency in the spatial autoregressive SF model.
On the other hand, Schmidt et al. (2009) were the first to make distributional assumptions on the two error terms that consider spatial dependence in the inefficiency error term. They specified the technical inefficiency error term using a prior distribution that captures unobserved spatial features and estimated the model through a Bayesian approach. Differently, Areal et al. (2012) introduced an autoregressive specification of the inefficiency error term to capture spatial dependence and used a Gibbs sampler and two Metropolis-Hastings steps to estimate the unknown parameters. Moreover, Fusco and Vidoli (2013) managed to solve the autoregressive error model proposed by Areal et al. (2012) using a maximum likelihood estimation technique instead of a Bayesian approach. Afterwards, Tsionas and Michaelides (2016), making use of a Bayesian estimator, developed a spatial SF model splitting the inefficiency error term into an idiosyncratic term and a component capturing regional spillover effects. Finally, Orea and Alvarez (2019) considered crosssectional dependence in both the inefficiency and the random error term and obtained a spatial SF model with closed form for the likelihood function using a transposed version of the Wang and Ho (2010) model. Specifically, Wang and Ho (2010) took advantage of the scaling property proposed by Wang and Schmidt (2002) to remove the time-invariant effects from the firm-specific inefficiency error term. Transposing the within temporal transformation to adapt to a cross-sectional setting, Orea and Alvarez (2019) replaced the firm-specific error term with a time-varying industry-specific term in order to develop a spatial SF model considering cross-sectional correlation in the inefficiency component that can be estimated using maximum likelihood techniques.
Considering the pros and cons of the two different modelling approaches, the most interesting advantage of introducing the spatial lag of Y and of X in SF models relates to the possibility of directly measuring direct and indirect effects originating from the variables that determine firms' productivity level and that also affect neighbouring producers. However, ignoring spatial dependence in the inefficiency error term can result in biased estimates of the inefficiency distribution and in different rankings of inefficiencies across agents (Areal et al., 2012;Schmidt et al., 2009). Moreover, estimating an SF model without taking cross-sectional dependence in the error term into consideration decreases the statistical efficiency of the model (Orea & Alvarez, 2019). On the other hand, considering spatial dependence in the error term is useful to take unobserved but spatially correlated variables such as environmental and climatic conditions or locationspecific attributes into account. Moreover, introducing a spatial structure in the inefficiency error term allows capturing spatial correlation deriving from the specific attributes that commonly characterize all the producers belonging to the same area and that affect firms' efficiency level. Nevertheless, not introducing the spatial lag of the dependent variable in SF models can result in biased estimates due to an omitted variable bias (Glass et al., 2016). Therefore, in this paper we merge all the different spatial structures described before in order to jointly consider all the different features outlined above, obtaining a comprehensive specification that takes four different sources of spatial dependence into account. Specifically, we introduce the spatial lag of Y and of X to capture global and local spatial spillovers affecting the frontier function as introduced by Glass et al. (2016), and we also add a spatial structure to the inefficiency and to the error term to capture behavioural and environmental spatial correlation, respectively, as proposed by Orea and Alvarez (2019). Hence, we obtain a spatial Durbin SF model for panel data introducing cross-sectional dependence both in the inefficiency and in the error term (SDF-CSD). The most appealing feature of our model is that it allows capturing global and local spatial spillover effects while controlling for spatial correlation related to firms' efficiency and to unobserved but spatially correlated variables. The SDF-CSD model can be estimated using maximum likelihood techniques, modifying the estimation procedure suggested by Orea and Alvarez (2019) in order to consider the endogeneity deriving from the inclusion of the spatial lag of the dependent variable. In the simulation section, we show that our spatial estimator is able to distinguish between frontier and error-based spillovers considering sparse spatial weight matrices (as binary contiguity or truncated inverse-distance matrices).
The remainder of the paper is structured as follows. The next section introduces the model specification, estimation procedure, and derivation of the marginal effects and of the technical efficiency (TE) scores. The third section tests the finite sample properties of our estimator by performing some Monte Carlo simulations. The fourth section shows a brief application to aggregated data from the Italian agricultural sector. The fifth section concludes.

Model specification
The spatial Durbin SF production function for panel data allowing for spatial cross-sectional dependence in both the inefficiency and the error term (SDF-CSD model) is specified as in equations (1-5), for i = 1, . . . , N and t = 1, . . . , T : where Y it represents the production output (usually expressed in log-form) of firm i at time t and X it is a 1 × k vector of production inputs used by firm i at time t with associated parameter vector b (k × 1). To take spatial dependence related to the frontier function into consideration, we include in the model the spatial lag of the dependent variable and the spatial lag of the production inputs, detecting global and local spatial spillovers through r (1 × 1) and u (k × 1), respectively. In particular, w ij represents the generic element of the spatial weight matrix W , built as a rownormalized (N × N ) matrix identifying the spatial structure of the cross-sectional units. Focusing on the two error terms, we follow the specification proposed by Orea and Alvarez (2019) for the definition of both the inefficiency and the random error component assuming that they are independently distributed across the cross-sectional dimension. In particular, the vector v t = (ṽ 1t , . . . ,ṽ Nt ) represents the error term distributed as a multivariate normal random variable with variance-covariance matrix P accounting for unobserved but spatially correlated variables thought the inclusion of M g , as defined in equation (3). Finally, assuming that the scaling property (Wang & Schmidt, 2002) holds, the inefficiency error term u it is decomposed as the product of two components: a scaling functionh it and a basic distribution u * t . Specifically,h it depends on a positive function f of firm exogenous variables Z it with associated parameter vector f as specified in equation (4), and on the spatial lag (I N − tW ) −1 capturing spatial dependence in the variables that determine technical efficiency. On the other hand, u * t is an industry-specific inefficiency error term following a truncated normal distribution with mean zero and variance s 2 u . Therefore, inefficiency depends on an industry-specific term common to all firms but varying in time, and on firm-specific exogenous variables that also affect neighbouring producers. Therefore, this model allows both the technical inefficiency error term u it =h it u * t and the random noiseṽ it to be cross-sectionally (spatially) correlated.
The use of the scaling property in defining the inefficiency error term is fundamental to obtaining a closed-form for the likelihood function. Indeed, considering cross-sectional dependence in the inefficiency error term generally precludes using standard maximum likelihood techniques for the estimation of the unknown parameters. Therefore, defining u it using the scaling property helps in overcoming this issue and allows to estimate the SDF-CSD using ML algorithms implemented in standard software.
For the model specified in equations (1-5), the following conditions must hold: (1) for a symmetric W , the spatial autoregressive parameters r, t and g must be within the range (1/v min , 1), where v min is the smallest eigenvalue of the spatial weight matrix W and the upper bound equals 1 if W is row-normalized; and (2) the cross-sectional correlation must converge to zero when the distance separating two spatial units goes to infinity (Kelejian & Prucha, 1998, 1999. The first condition ensures that (I N − rW ), (I N − tW ) and (I N − gW ) are non-singular matrices, and the second one guarantees that the row and columns sums of W , (I N − rW ) −1 , (I N − tW ) −1 and (I N − gW ) −1 , before W is row-normalized, are uniformly bounded in absolute value as N tends to infinity. In particular, assumption (2) is always satisfied when W is a binary contiguity matrix while for inverse-distance matrices it should be introduced a cut-off point in the distance measure (Elhorst, 2010).

Estimation procedure
The partial loglikelihood function, referring to a single time period t, is shown in equation (6), where m * and s 2 * are defined as in equations (7) and (8),h t = (h 1t , . . . ,h Nt ), and the generic element 1 it belonging to the vector 1 t is specified as in equation (9). Moreover, F represents the cumulative distribution function of the standard normal random variable and log |I N − rW | is the determinant of the Jacobian of the transformation from 1 it toY it . Indeed, since ∂1 ∂Y = (I N − rW ), we take into account the endogeneity deriving from the inclusion of the spatial lag of the dependent variable. Starting from equation (6), the final loglikelihood function is equal to L(Q; y) = T t=1 ℓ t , that is, the sum of the partial loglikelihood functions computed for each year in the sample: Consistent parameter estimates can be obtained by maximizing the final loglikelihood function using a numerical constrained maximization algorithm. 1 According to the first assumption described above, the parameter space for an autoregressive process is (1/v min , 1), where v min is the smallest eigenvalue of W and the upper bound equals 1 because W is row-normalized (Kelejian & Prucha, 1999). Therefore, all the autoregressive parameters (r, t, g) should be bounded to the previous interval to ensure that (I N − rW ), (I N − tW ) and (I N − gW ) are non-singular matrices. Moreover, the two variance parameters s 2 u and s 2 v should only take positive values. A further characteristic of the SDF-CSD model is that its loglikelihood function is the same both considering a production function and a cost function. For more details on the estimation procedure, see Appendix A in the online supplementary materials.
Following Orea and Alvarez (2019), technical efficiency scores (TE) can be found substituting the parameter estimates in TE it = exp ( − E(u it |1 it )), where the conditional expectation of u it given1 it is shown in equation (10):

Marginal effects
As usual in spatial models including the spatial lag of Y , the b estimates cannot be interpreted as marginal effects because changes in the generic regressor X r of firm i influence the production output of firm j. To show this, the SDF-CSD model is written in the reduced form in equation (11), where Y is an (NT × 1) vector representing firms' output, X is an (NT × k) matrix containing the k production inputs with associated parameter vector b (k × 1); W is the block diagonal NT × NT spatial weight matrix; I NT is the NT × NT identity matrix;ṽ is the error term distributed as shown in equations (2) and (3); and u is the inefficiency error term specified following the scaling property as described above. Moreover, r is the autoregressive parameter capturing global spatial spillovers; u is the parameter vector (k × 1) associated with input spillovers; t, belonging to the spatial structure included in the scaling function, is the parameter capturing behavioural correlation depending on spillovers effects related to the determinants of firms' inefficiency; and g, embedded in the structure ofṽ, is the spatial parameter capturing environmental spatial dependence associated to unobserved but spatially correlated variables.
Deriving equation (11) with respect to the generic regressor X r leads to Thus, it is evident that in spatial models that include that spatial lag of Y , the b estimates no longer represent the marginal effects of the X variables on Y . LeSage and Pace (2009) proposed to compute the marginal effects starting from the S r (W ) matrix differentiating into (1) direct effects computed as the average of the diagonal elements of S r (W ); (2) indirect effects calculated as the average of the sum of the non-diagonal elements of S r (W ); and (3) total effects equal to the sum of the direct and the indirect effects. The related standard errors or t-statistics can be computed either following the simulation method proposed by LeSage and Pace or using the delta method.
In order to compute the marginal effect of a generic Z variable on firms' inefficiency level, we consider the unconditional definition of firms' inefficiency. Indeed, within our SDF-CSD framework, the expression t provides a measure of u it that is conditional on the frontier feedback effects. Starting from the reduced form in equation (11), the unconditional definition of the inefficiency error term is given byũ it = (I N − rW ) −1 u it . Thus, the first derivative ofũ it with respect to Z, representing the marginal effect of Z onũ it , includes the spatial filter (I N − rW ) −1 related to the endogenous spatial lag of Y , as shown in equation (12a-d). However, it should be noted that both the conditional and the unconditional definitions of the inefficiency error term lead to the same interpretation of f.
In particular, we specify f (Z it , f) = exp (Z it f) and we apply the logarithm toũ it , as suggested by Wang and Schmidt (2002), to generate a very simple expression for the effect of the generic determinant Z on firms' inefficiency level. Thus, the interpretation of f is straightforward, as it represents the marginal effect of Z on log (ũ it ), that is, the logarithm of the level of technical inefficiency of the firm.

MONTE CARLO SIMULATIONS
To test the finite sample properties of our spatial estimator, we simulated NT data. In particular, we defined for each time period t: the input variable X and the exogenous variable Z that determines technical inefficiency as standard normal random vectors (N × 1), h as an exponential function and W as an (N × N ) row-normalized spatial weight matrix. Moreover, to introduce cross-sectional correlation in the error term, we defined v as a multivariate normal random vector (N × 1) with zero mean and variance-covariance matrix equal to P, with Finally, the inefficiency error term u * , common to all firms but varying in time, is generated as a random value drawn from a truncated normal distribution with zero mean and variance equal to s 2 u . Therefore, following the scaling property, the complete inefficiency term u (N × 1) accounting for cross-sectional dependence, is equal to the product of u * , of the positive function of the firm exogenous variable h = exp (Zf), and of the spatial filter (I N − tW ) −1 . Then, the dependent variable Y , for each time period t, is generated as a (N × 1) vector, as defined in equation (13), and the whole procedure is repeated for all time periods. Lastly, the partial loglikelihood function ℓ t is calculated for each time period and the final loglikelihood function is obtained as We performed three Monte Carlo simulations to test how the choice of different spatial weight matrices, of different true values of the parameters and of the sample sizes affects the bias, the standard deviation (SD) and the mean squared error (MSE) of the estimated parameters. In particular, in the first simulation, we let N and T vary across different values (N = 100, 200, 300 and T = 5, 10, 15), keeping the true values of the parameters constant, while in the second one we change the true values of the parameters keeping N and T fixed at N = 100 and T = 10. Moreover, in the third block of simulations, we consider different spatial weight matrices. In particular, in the first two blocks of simulations, W is defined as a binary contiguity spatial weight matrix, while in the third block of simulations, also inverse-distance spatial weight matrices with different truncation criteria are taken into consideration, keeping the true values of the parameters and N and T constant. In the first and third blocks of simulations the true values of the parameters are fixed at {b = 0.50, r = 0.30, u = 0.30, f = 0.50, t = 0.30, g = 0.30, s 2 u = 0.10, s 2 v = 0.20}. All the simulations are performed using 1000 repetitions. Table 1 shows the results of the Monte Carlo experiment for different values of N and T , Table B1 in Appendix B in the online supplementary materials shows how the final sample properties of the estimated parameters are affected by changes in one value at a time, keeping N and T constant, while Table 2 takes different kinds of spatial weight matrices into consideration with T and N fixed at 5 and 300, respectively.
The results in Table 1 show that b, u, g, s 2 u and s 2 v are estimated with negligible bias even using a very small sample while the bias of the other parameters (r, f, t) quickly reduces increasing N or considering more time periods. In particular, r, f and t tend to be slightly underestimated with very small N and T but overall, the bias of all parameters approaches zero considering a sufficiently high sample size. Likewise, the SD and the MSE tend to decrease, increasing the numerosity of the sample. Moreover, Table B1 in Appendix B in the online supplementary materials, presenting the results of the second block of simulations for T = 10 and N = 300, shows that the estimates are robust to different changes in the true value of one parameter at a time, keeping all the other constant.
Finally, Table 2 shows the results of the simulations considering different kinds of spatial weight matrices. Differently from the previous simulations taking a binary contiguity spatial weight matrix into consideration, here, W indicates a dense inverse-distance spatial weight matrix, W50 and W30 indicate two inverse-distance spatial weight matrices truncated at 50 and 30 km, respectively, and W250n, W100n and W50n stand for three inverse-distance spatial weight matrices considering only the 250, 100 and 50 nearest neighbours, respectively. All these spatial weight matrices have been created starting from a sample of 300 tourism facilities located in the North-West of Italy belonging to the ATECO55 sector. Moreover, all these matrices are row standardized, T is fixed at 5, and the number of replications is equal to 1000.
The results shown in Table 2 indicate that the bias of b, r, u, f and of the two variance parameters s 2 u and s 2 v is not affected by the choice of W. Conversely, t and g tend to be underestimated if a large number of neighbours is considered. Indeed, the bias of t and g is equal to −0.0661 and −0.0206 choosing a dense W, while it decreases to −0.0172 and −0.0079 truncating W at 30 km. Moreover, also the standard deviations associated with these two parameters tend to be quite high considering a dense W but they tend to decrease for W50 and W30. Similarly, t and g tend to be underestimated when W is defined as an inverse-distance spatial weight matrix considering the n nearest neighbours and the number of neighbours is high, but the bias decreases as n reduces. Indeed, using W250n the bias associated with t and g is equal to −0.0545 and −0.0131 respectively, but it reduces to −0.0283 and −0.0034 for W50n. Therefore, to obtain unbiased and more efficient estimates for the unknown parameters it is better to consider sparse matrices such as a binary contiguity spatial weight matrix or an inversedistance matrix with a truncation point. This issue has already been addressed by Mizruchi and Neuman (2008) that suggested preferring sparse spatial weight matrices when using spatial models. Indeed, in accordance with our simulation results, they demonstrated that dense inverse-distance matrices are likely to produce negatively biased parameter estimates and that this bias becomes more severe at higher levels of network density.
Alternatively, if a dense W is anyhow preferred, to cut down the bias detected in this simulation study, it can be assumed that the spatial weight matrix is not identical among the different spatial structures considered in the SDF-CSD model. Therefore, we implement a further simulation considering two different spatial weight matrices to evaluate if the bias detected using a dense W diminishes. The results, shown in Table B2 in Appendix B in the online supplementary materials, indicate that using different spatial weight matrices to differentiate among the different kinds of spatial effects is sufficient to cut off the bias of g, while to obtain an unbiased estimate of t it is necessary to consider a limited number of neighbours. In conclusion, the only limitation in estimating the SDF-CSD model without bias is to necessarily assume that the spillover effects associated with the inefficiency error term result from closest neighbours. This assumption should not be considered too restrictive because spatial dependence influencing neighbouring firms' efficiency level mainly arises from local factors such as emulation, face-to-face interactions, local cooperation, and individual contact (Griliches, 1992).  We take advantage of the SDF-CSD model to analyse the performance of the Italian agricultural sector considering four different kinds of spatial effects. This novel spatial estimator results to be particularly suitable for this sector due to the strong importance of unobserved location-specific attributes in the agricultural industry. Indeed, thanks to the spatial structure attached to the random error term it is possible to consider spatial cross-sectional dependence arising from unobserved factors common to neighbouring areas such as soil conditions or climatic, topographic and environmental characteristics (Schmidt et al., 2009). Moreover, the spatial structure related to the determinants of firms' efficiency allows capturing behavioural spatial dependence arising from emulation behaviours of agricultural producers located in nearby areas and from policies and institutions operating at the local level (Areal et al., 2012). On the other hand, the spatial structures entering in the frontier function relate to productivity and input spillovers. Input spillovers may generate from a greater availability of specific products, input suppliers, assets and workers with industry-specific skills in a certain area (Marshall, 1890). Furthermore, farmers' productive performance may be related to the one of neighbours due to the transmission of knowledge and best practices between peers (Cardamone, 2020), collective behaviours such as similar financial decisions resulting from faceto-face relationships, exchange of ideas and learning from others (Skevas & Lansink, 2020), farmers' adoption of new similar technologies to face specific techno-economic problems that are common to firms operating in nearby areas (Billé et al., 2018), and marketing-related externalities such as positive feedbacks deriving from 'protected designation of origin' (PDO) certifications (Vidoli et al., 2016). Indeed, the successful performance of neighbouring producers may generate economic returns also to peers because of the increased reputation of the whole area (i.e., 'halo effect', Beebe et al., 2013). Hence, in this section, we estimate the SDF-CSD model using data on the Italian agricultural sector at NUTS-3 level in the time period 2008-18. The data are collected from the Rete di Informazione Contabile Agricola (RICA) survey, which is the Italian counterpart of the Farm Accountancy Data Network (FADN) survey. The FADN survey represents the only harmonized source of microeconomic data on the evolution of incomes and on the economic-structural dynamics of farms at the European level. The FADN sample does not represent the entire universe of farms surveyed in a given territory, but only those that, due to their economic dimension 2 can be considered professional and market oriented. The methodology adopted aims to provide representative data on three dimensions: region, economic dimension, and technical-economic order. In this empirical application, we consider aggregated data at NUTS-3 level from the RICA survey because for confidentiality reasons it is not possible to know the exact location of each farm in the sample. Moreover, aggregated data at the municipal level do not ensure representative results because in most cases they do not have the minimum required sample numerosity of 5 units. Thus, the NUTS-3 aggregation results to be the finer aggregation level that also guarantees representative results.
We use a Cobb-Douglas specification to model the production function equation, as shown in equation (14) for i, j = 1, . . . , N (i = j) and t = 1, . . . , T . Since we include four input variables in the analysis, we prefer a Cobb-Douglas functional form with respect to a Translog specification because it involves the estimation of fewer parameters and thus, it facilitates the interpretation of the results. Moreover, the Cobb-Douglas function is often used to estimate the production function parameters due to its ability to provide more efficient estimates, especially when dealing with small samples (Yao & Liu, 1998).
Specifically, the dependent variable Y it is defined as the logarithm of the total value added generated by the agricultural sector in province i at time t. Following Vidoli et al. (2016), we consider four input variables, all in log-form: total working hours (L it ), utilized agricultural area (AA it ), machinery (M it ), and water, energy and fuel (WEF it ). Moreover, we add the time trend t to the model specification to capture the temporal dynamics, where t has a minimum value of 1 for 2008 and a maximum of 11 for 2018. To consider cross-sectional spatial dependence affecting the frontier function we include in equation (14) the spatial lag of the dependent variable and of the four production inputs. Specifically, following the queen contiguity criteria, we identify neighbouring observations as provinces sharing a common border or vertex. Thus, we build the spatial weight matrix W as a row standardized binary contiguity matrix considering second-order neighbours (i.e., we assign value 1 to neighbouring provinces and to neighbours of neighbours and 0 otherwise 3 ). In particular, w ij refers to a generic spatial weight belonging to the (N × N ) sparse spatial weight matrix having all zeros on the main diagonal. Hence, we are able to detect global spatial dependence through r, while the parameters u L , u AA , u M , and u WEF measure local spatial dependence related to the four input variables. Finally,ṽ it and u it represent the two error terms. The random disturbance accounts for unobserved but spatially correlated variables thought g, capturing remaining spatial dependence resulting from unobserved but cross-sectionally correlated spatial features such as climatic conditions, soil characteristics, and local institutional and socio-economic factors. On the other hand, the inefficiency error term, defined using the scaling property, can be written as u it =h it u * t , whereh it is the scaling function defined as in equation (15) and u * t is the industry-specific inefficiency error term common to all firms but varying in time following a truncated normal distribution with mean 0 and variance s 2 u .
The scaling function shown in equation (15) is composed by the spatial lag (I N − tW ) −1 and by a positive function of the determinants of farms' inefficiency level. As usual in SF model specifying the inefficiency error term using the scaling function approach, we choose the exponential function to obtain easily interpretable estimates for the f parameters. Moreover, as proposed by Orea and Alvarez (2019), we are able to capture the overall spatial dependence associated with the variables that determine cross-sectional inefficiency through t. Indeed, farms' inefficiency level can also be related to the observable characteristics of neighbouring producers and ignoring these features can lead to heteroskedasticity issues in the inefficiency error term.
Focusing on the determinants of inefficiency, we model the scaling function introducing some variables related to the main characteristics of Italian farms at the provincial level. Specifically, we consider the average farm dimension in the province introducing two variables (Small and Big) that respectively measure the proportion of small and big farms. We take the organizational type of the farms into consideration including in the model the variable Family that equals the proportion of farms run by family members in the province. Moreover, we include in the scaling function equation the variable Diversified measuring the proportion of diversified farms in the province and Hired measuring the average percentage of hired land. Furthermore, farmers' characteristics such as age, gender, and education are often considered as a proxy for management practices and social capital. Therefore, we include in the scaling function Youth and Woman representing the percentage of farms run by young entrepreneurs (with less than 40 years) and by women, respectively. Finally, we investigate if subsidies positively or negatively affect farmers' inefficiency levels. Thus, we introduce the variable Subsidies defined as the amount of total subsidies received by farms located in each province over total farm income as proposed by Stetter and Sauer (2021). Further details on the definition and measurement of the variables considered in the analysis and some brief descriptive statistics can be found in Table 3. Table 4 shows the estimation results of the SDF-CSD model (yxuv model), of the SF models introducing only frontier and error-based spillovers (yx and uv model, respectively) and of the non-spatial SF model. The complete estimates of all the nested specifications are contained in Table C1 of Appendix C in the online supplementary materials. Indeed, starting from the SDF-CSD it is possible to restrict the model specification in order to capture only specific sources of spatial dependence depending on the economic phenomenon under investigation. As a consequence, the SDF-CSD model leads the way to a number of spatial specifications never introduced before but that can be very useful in empirical applications. Moreover, it is worth highlighting that we did not have any problem simultaneously estimating the four different spatial parameters using a numerical maximization algorithm even if our sample numerosity is not very large. Comparing the estimation results of the nested models shown in Table 4, it can be noticed that the b estimates and the f estimates are quite robust to the different specifications. Nevertheless, as previously described, the b coefficients cannot be interpreted in a meaningful way when the spatial lag of Y is included in the model (a brief description of the marginal effects of the input variables is presented in subsection 4.3). As for the f estimates, our results indicate that provinces with a higher percentage of small farms as well as provinces with a lot of farms run by family members tend to be more inefficient while big farms contribute to decreasing the inefficiency level of the agricultural sector at NUTS-3 level. Moreover, we find that the degree of diversification positively affects inefficiency. Indeed, specialization might boost technical efficiency more than diversification since it enables farmers to specialize in a few tasks, and therefore it improves management practices (Latruffe, 2010). Investigating the link between technical inefficiency and the use of external factors, we find a negative effect of the percentage of rented land in the province on technical inefficiency even if Hired results to be significant only in some  nested specifications. Considering subsidies, our results confirm the idea that they contribute to decreasing the efficiency level of farms likely due to lowered farmers' motivation and distorted farms' production structure and factor use (Rizov et al., 2013). Finally, provinces with a higher percentage of farms run by women or by young entrepreneurs tend to be more inefficient than others. Indeed, as Chen et al. (2009) and Mathijs and Vranken (2001) explained, older farmers result to be more efficient thanks to increased experience and learning ability. Concentrating on the spatial parameters, we detect positive and significant spatial dependence either at the frontier level or related to the two error terms. Thus, the estimate of r indicates that positive global productivity spillovers affect the Italian agricultural sector as well as positive spillover effects related to the determinants of farms' inefficiency (t), and to unobserved environmental factors (g). Specifically, the degree of global spatial dependence (0.14) estimated using the SDF-CSD model is considerably smaller than the level of behavioural and environmental spatial correlation associated with the two error components (both equal to 0.37) indicating that the Italian agricultural sector is more strongly affected by error-based than by frontierbased spatial dependence.

Results
Finally, in Table C2 of Appendix C in the online supplementary materials we compare the different nested models using some likelihood ratio tests and information criteria. According to both the LR test and to the AIC criteria, in this case study, it is better to simplify the model specification using a Manski model that, besides considering the spatial lag of Y and of X , includes one unique spatial term referred to the whole error structure. Indeed, the degree of spatial dependence associated with the inefficiency and the random error term is the same for both the two error components using the SDF-CSD model.

Robustness check and marginal effects
As a robustness check, we estimate the SDF-CSD model using different spatial weight matrices. Besides the second-order binary contiguity matrix (W2), we consider a first-order binary contiguity matrix (W) and a dense inverse-distance matrix (Wid). According to the results of the Monte Carlo simulations, when we introduce the dense inverse-distance matrix we associate it only to the spatial lags related to the frontier function while for the two error terms we use a sparse matrix like W or W2. The results of the robustness check are shown in Table E1 of Appendix E in the online supplementary materials while Table 5 shows the corresponding marginal effects.
The direct effects of the X variables on Y do not vary depending on the choice of the spatial weight matrix. Specifically, they are all positive and significant, and labour results to be the more effective input variable (0.61), followed by water, energy, and fuel (0.20). In line with the results of Fusco and Vidoli (2013), capital inputs such as machinery (0.07) and land (0.13) contribute less to the productive performance of the Italian agricultural sector compared with human resources. Overall, the sum of the elasticities related to the four input variables results to be greater than 1 regardless of the spatial structure considered. Thus, in line with most of the European countries (Rizov et al., 2013), the Italian agricultural sector is characterized by increasing returns to scale, indicating that increasing the inputs by 1% would produce more than a 1% increase in output. Considering the indirect effects, they tend to vary as the choice of the spatial weight matrix changes. In particular, the indirect effect associated with labour is positive and significant only when we consider a first-order contiguity matrix related to the frontier function. Conversely, the indirect effect of machinery is positive and significant only considering a dense inverse-distance spatial weight matrix associated with the frontier function. Conversely, the indirect effect of utilized agricultural areas and water, energy, and fuel does not change across the different estimated models. Indeed, while the indirect effect of AA is always negative and significant, the indirect effect of WEF is never significantly different from zero. In sum, we detect positive spillover effects at the local level in terms of labour while the positive spillover effects Note1: * * * :pvalue ≤ 0.01; * * :pvalue ≤ 0.05; * :pvalue ≤ 0.10. Note2: W2, second-order contiguity matrix associated with all the spatial terms; W, first-order contiguity matrix associated with all the spatial terms; W2,W, second-order contiguity matrix associated with the frontier function and first-order contiguity matrix associated with the two error components; W,W2, first-order contiguity matrix associated with the frontier function and second-order contiguity matrix associated with the two error components; Wid,W2, dense inverse-distance matrix associated with the frontier function and second-order contiguity matrix associated with the two error components; and Wid,W, dense inverse-distance matrix associated with the frontier function and first-order binary contiguity matrix associated with the two error components.
A spatial stochastic frontier model including both frontier and error-based spatial cross-sectional dependence related to machinery result to be significant only at the global level. On the contrary, provinces that dispose of larger agricultural areas tend to be located near provinces with fewer lands likely due to the specific territorial conformation and the alternation of more rural provinces with more urbanized ones.

4.4.
Technical efficiency scores Figure 1 shows the kernel density plots of the technical efficiency scores computed estimating the SDF-CSD model, the non-spatial SF model, and the uv and xy models, introducing spatial dependence only in the error terms and in the frontier function, respectively. Considering the scale of the distribution, it can be noticed that the TE scores' distribution from the non-spatial model is the closest to the one obtained through the SDF-CSD. On the other hand, the xy model tends to overestimate the TE scores while the uv model tends to underestimate them. As for the shape of the distribution, the distribution obtained using the SDF-CSD model resembles the one coming from the spatial error model but is mitigated by the quasi-normal shape of the xy model. Thus, the SDF-CSD model allows to accurately combine the main features of the distributions of the TE scores resulting from the uv and xy model. In conclusion, this empirical application shows that the non-spatial TE scores distribution is the one that better approaches the SDF-CSD one while considering spatial dependence only in the frontier function or in the error terms leads to severe distortions both in the scale and in the shape of the distribution. Therefore, it is fundamental to combine frontier based and error-based spatial effects both to consider different sources of spatial dependence through a comprehensive model and also to obtain consistent estimates of the TE scores. Moreover, in Table D1 of Appendix D in the online supplementary materials we present the efficiency ranking of the 107 Italian provinces where it can be observed more in detail how both the ranking and the level of the TE scores modify depending on the kind of spatial structures considered. In addition, in Table D2 of Appendix D we show how the efficiency score ranking changes in time by comparing the 2008 and 2018 results. Overall, it can be observed that the  Figure 2 maps the TE scores resulting from the SDF-CSD model for the 107 Italian provinces, comparing 2008 and 2018. Overall, provinces located in the North of Italy result to be much more efficient than those located in the South and this gap has remarkably increased in the past eleven years. Specifically, the large efficiency cluster in the Po Valley is strengthening over time while Sardinia, Calabria, the southern Apulia, and the Northern provinces at the border with Switzerland are achieving lower and lower efficiency scores. Besides highlighting a severe North-South divide that is increasing with time, Figure 2 shows that the Italian provinces are characterized by strong spatial concentration considering the efficiency level of the agricultural sector.
Further insights on the spatial features related to technical efficiency scores are contained in Appendix D in the online supplementary materials.

CONCLUSIONS
In this paper we propose a novel spatial SF model that introduces four different kinds of spatial spillover effects, two related to the frontier function (i.e., productivity and input spillovers) and two related to the error structure (behavioural correlation associated with the inefficiency error term and environmental correlation referring to unobserved variables). While previous literature only considers frontier based or error-based spillover effects, this is the first work that merges together the two different approaches, obtaining a very comprehensive specification never introduced before. Despite his complex spatial structure, one of the most remarkable features of the SDF-CSD model is that, thanks to the modelling approach suggested by Orea and Alvarez (2019) for the two error terms, it can be estimated using standard maximum likelihood algorithms and it allows to straightforwardly interpret the effects of the inefficiency determinants. Investigating whether our model is able to distinguish in practice between the four different sources of spatial dependence, we showed that our spatial estimator is able to differentiate between frontier and error-based spillovers considering sparse spatial weight matrices (as binary contiguity or truncated inverse-distance matrices). Alternatively, unbiased estimates can be obtained using two different spatial weight matrices simultaneously and associating a sparse matrix to the spatial lag of the inefficiency error term. Therefore, the only limitation in estimating the SDF-CSD model without bias is to necessarily assume that the spillover effects associated with the inefficiency error term result from closest neighbours. This assumption should not be considered too restrictive because geographical proximity is fundamental for interaction, cooperation and for the transmission of new knowledge (Griliches, 1992).
In the empirical section of this paper, we applied the SDF-CSD model to NUTS-3 level data on the Italian agricultural sector. The main advantage resulting from estimating the SDF-CSD model in practice concerns the possibility of also evaluating a number of nested specifications never introduced before. Indeed, starting from the SDF-CSD and making some LR tests for nested models, it is possible to test whether it is better to simplify the model specification considering only specific spatial lags or if a comprehensive spatial SF model is required. Thus, it is possible to precisely assess which kind of spatial effect is more appropriate for studying the phenomenon under investigation without making a priori assumptions on the spatial structure of the data. Finally, we showed that the SDF-CSD model allows obtaining more accurate technical efficiency scores, mixing together the main distributional features resulting from spatial SF models that only include frontier based or error-based spillovers.
From a practical perspective, the empirical evidence on the existence of significant agglomeration externalities originating from various channels has important implications for policy makers dealing with the Italian agricultural sector. Indeed, strategic plans and programmes aimed at improving the Italian agricultural sector's performance may take advantage of positive spillovers by strengthening the cohesiveness of the networks, providing more opportunities for farmers to learn from each other, stimulating cooperation and networking, and encouraging local farmers' associations as well as knowledge dissemination on the use of new equipment and new management practices. Exploiting existing spatial interactions by creating a collaborative environment can be an effective strategy for policy makers to overcome the technological backwardness of most Italian farmers and thus, to boost the productivity of the Italian agricultural sector.
Given the completeness of the SDF-CSD model in considering the different spatial structures, in future extensions of this work, it could be interesting to modify or relax the underlying modelling assumptions. In particular, it could be useful to introduce the possibility to account for cross-sectional variation in the inefficiency error term by implementing a Bayesian algorithm for the estimation. Moreover, rather than assuming that the inefficiency distribution is half-normal, different specifications for the two error components can be considered, such as an exponential or a gamma distribution.

DISCLOSURE STATEMENT
No potential conflict of interest was reported by the author.

NOTES
1 Matlab codes are available from the author upon request. 2 In Italy, the minimum threshold for the inclusion in the RICA observation sample corresponds to a minimum standard gross income of €8000. 3 A second-order matrix compared with a first-order one allows one to consider more complex spatial structures and better identify spatial clusters. Moreover, this matrix minimizes the loglikelihood function. Some robustness checks using different spatial weight matrices are provided afterwards.