Inference in the spatial autoregressive efficiency model with an application to Dutch dairy farms

This article extends the conventional spatial autoregressive efficiency model by including firm characteristics that may impact efficiency. This extension allows performing the typical inference in spatial autoregressive models that involves the derivation of direct and indirect marginal effects, with the latter revealing the nature and magnitude of spatial spillovers. Furthermore, this study accounts for the endogeneity of the spatial autoregressive efficiency model using a lag spatial lag efficiency component, which makes inference to be performed in a long-run framework. The case study concerns specialized Dutch dairy farms observed over the period 2009–2016 and for which exact geographical coordinates of latitude and longitude are available. The results reveal that the efficiency scores are spatially dependent. The derived marginal effects further suggest that farmers’ long-run efficiency is driven by changes in both their own and their neighbors’ characteristics, highlighting the existence of motivation and learning domino effects between neighboring producers. © 2019 The Author. Published by Elsevier B.V. This is an open access article under the CC BY license. ( http://creativecommons.org/licenses/by/4.0/ )


Introduction
A plethora of studies in the economics/operational research literature is engaged with the measurement of producers' inefficiency, thus relaxing the behavioral assumption that they are perfect decision-making units. In a parametric setting, the technique of Stochastic Frontier Analysis (SFA) introduced by Aigner, Lovell, and Schmidt (1977) and Meeusen and van den Broeck (1977) is almost exclusively used to measure producers' performance. The main advantage of SFA lies on that deviations from optimality are not only attributed to pure inefficiency but also to factors that are outside the control of producers. Measuring producers' performance and identifying its potential drivers can provide important information to both businesses and policy-makers. Firms can use this information to improve their efficiency levels, while policymakers can get informed about the success of their instruments in improving producers' performance.
Despite the numerous advancements in the measurement of firms' efficiency that allow drawing such important conclusions for both producers and policy-makers, an important aspect that has largely been ignored in the related literature is the potential existence of spatial dependence among producers. The first law of E-mail address: ioannis.skevas@ucc.ie geography introduced by Tobler (1970) states that "everything is related to everything else, but near things are more related than distant things". Bockstael (1996) and Weiss (1996) discussed this issue in the context of agriculture arguing that economic processes such as agricultural production are spatial phenomena. This may be highly relevant when measuring producers' performance as spatial dependence in their efficiency scores may be observed. For instance, farmers' motivation may be influenced by that of their neighbors and result in changes in their working effort and/or investment choices and consequently their efficiency levels ( Areal, Balcombe, & Tiffin, 2012 ). Furthermore, if neighboring farmers communicate with each other and exchange information regarding their production practices, individuals may learn how to use their resources more efficiently by collaborating with more experienced neighbors and therefore improve their efficiency ( Pede, Areal, Singbo, McKinley, & Kajisa, 2018 ).
The hypothesis regarding the existence of spatial dependence in producers' efficiency levels has only recently been tested by few empirical studies. The aim of early studies on spatial dependence of producers' efficiency levels was to identify whether efficiency scores follow a certain pattern dependent on their location. Examples of such studies include Hadley (2006) and Zhu, Karagiannis, and Oude Lansink (2011) who did not specify a standard spatial econometric model but simply included region-specific dummies as potential determinants of efficiency. Other studies such JID: EOR [m5G;November 14, 2019;9:28 ] as Tsionas and Michaelides (2016) and Carvalho (2018) used regional information to define neighboring producers and specified a spatial autoregressive efficiency model to test for spatial dependence in their efficiency levels. This model assumes that individuals' efficiency depends on their neighbors' efficiency and a noise component 1 . The only studies that combined the use of the spatial autoregressive efficiency model with more micro-level geographical information are those of Areal et al. (2012) , Fusco and Vidoli (2013) and Pede et al. (2018) . The first two studies identified neighboring producers based on a given kilometre grid-square level and the municipality centroid, respectively, while the last study was the only one that had access to exact geographical coordinates (i.e. latitude and longitude). All the above studies reported a positive dependence in neighboring producers' efficiency scores. However, the present study argues that irrespective of the scale considered, the spatial autoregressive efficiency model employed by the above studies raises some theoretical as well as methodological challenges. Theoretically-wise, the aforementioned studies assume that individuals' efficiency is only related to their neighbors' efficiency and an error component, which is somewhat restrictive as one would expect that it would also be influenced by individuals' characteristics (i.e. subsidies received, management practices etc.). More importantly, neglecting these characteristics raises a methodological issue as it does not allow the above studies to perform the required inference and make statements regarding the nature and magnitude of spatial spillovers. According to Anselin (1988) , LeSage andPace (2009) andElhorst (2014) , inference in a spatial autoregressive model is not directly based on the parameter associated with the neighbors' dependent variable (in this case neighbors' efficiency), but on the derivation of the so-called direct and indirect marginal effects of the utilized firm characteristics. The direct effect measures the change in individuals' technical efficiency as a result of a change in one of their own characteristics, while the indirect effect quantifies the change in individuals' technical efficiency due to a change in one of their neighbors' characteristics. Therefore, the derivation of indirect effects allows one to make statements regarding the type and magnitude of spatial spillovers. For instance, if an increase in neighboring producers' experience contributes to an increase in individuals' efficiency, this implies that neighboring producers may indeed communicate and share their knowledge and experience, which ultimately leads to the individuals' efficiency gains. Finally, all the aforementioned spatial efficiency studies disregarded the presence of endogeneity in the spatial autoregressive model, which may result in biased estimates ( LeSage & Pace, 2009 ). Endogeneity exists due to simultaneity as the dependent variable of any individual appears both on left and the right-hand side of the spatial autoregressive model.

ARTICLE IN PRESS
This study adds to the literature that explores the existence of spatial dependence in producers' efficiency. Although it is interesting to merely study whether producers' efficiency levels are spatially dependent, it is more informative to know the nature and magnitude of spatial spillovers. If, for instance, empirical evidence shows that spatial spillovers exist because producers communicate with each other, policy-makers can take advantage of this information domino effect and promote more efficiently new technologies that can improve producers' performance. Specifically, this study contributes to the literature in four ways. First, it uses a farm-level panel dataset that contains information on the exact geographical location of farms based on their latitude and longitude. As mentioned above, studies that test for the existence of spatial dependence in producers' efficiency scores are almost exclusively based on regional information, thus not being able to differentiate between closer and more distant neighbors. Second, this study uses the typical spatial autoregressive efficiency model but further assumes that individuals' efficiency not only depends on neighbors' efficiency but also on certain characteristics. This is something that has totally been ignored by previous spatial efficiency studies. Third, the above addition also leads to a methodological contribution. Accounting for farm characteristics in the spatial autoregressive efficiency model allows this study to perform inference regarding the type of spatial spillovers and their magnitude. Fourth, this study tackles the issue of endogeneity in this spatial autoregressive modelling setting by using neighbors' lag efficiency levels, which makes inference being drawn in a long-run framework.
The proposed approach is applied to a very recent panel dataset of specialized Dutch dairy farms for which latitude and longitude information are available. This specific case study may be highly relevant for studying spatial spillovers among farmers since Groeneveld, Wesseler, and Berentsen (2013) argue that Dutch dairy farmers have the tendency to communicate with their neighbors and be aware of their production practices. The next section presents the modelling approach and the estimation method. A description of the data and the empirical specification follows. The results are then presented and the final section concludes.

Spatial stochastic frontier model
The stochastic frontier model has been systematically used to study the efficiency of agricultural decision-making units ( Emvalomatis, Stefanou, & Oude Lansink, 2011;Heshmati, Kumbhakar, & Hjalmarsson, 1995;Paul & Shankar, 2018;Reinhard, Lovell, & Thijssen, 20 0 0;Skevas, Emvalomatis, & Brümmer, 2018b ). The production frontier utilized in this study describes the maximum output that can be produced given the amount of inputs used ( Kumbhakar & Lovell, 2003 ). The model stacked over N individuals in period t is written as: where y t is a N × 1 vector of the logarithm of output in time t , a is a N × 1 vector of farm effects assumed to be draws from a Normal distribution with 0 mean and variance σ 2 α I , X t is a N × K matrix of the logarithm of K inputs in period t , β is the associated K × 1 vector of parameters to be estimated, v t is a N × 1 vector of two-sided error terms that account for statistical noise in time t and are assumed to be Normally distributed with 0 mean and variance σ 2 v I and TE t is the N × 1 vector of efficiency components in t period 2 . The efficiency term is projected from the unit interval to the real line (i.e. −∞ to + ∞ ) in order to avoid criticism that can be raised by subsequently imposing a lag spatial autoregressive process with environmental variables on a bounded variable such as efficiency (see Tsionas, 2006 ). The inverse of the logistic function is used for this purpose and the transformed technical efficiency is defined as s t = log ( TE t 1 −TE t ) . In a similar setting, the same transformation for the efficiency component was used by Emvalomatis (2012) who specified a autoregressive process on efficiency and Paul and Shankar (2018) who used an efficiency 2 Note that the stochastic frontier model presented in Eq.
(1) is parameterized in terms of the efficiency component TE t and not in terms of the typical ineffi- . This is because the subsequent need to transform this component so that it takes any values between −∞ and + ∞ can better be fulfilled for variables that lie on the unit interval than those that can take any non-negative value. This is because transformation of the former can be materialized with functions that are more prominent in the econometrics literature (see probit and logit models). JID: EOR [m5G;November 14, 2019;9:28 ] effects model. Subsequently, one can follow previous spatial efficiency studies and assume that individuals' transformed efficiency depends on neighbors' transformed efficiency and a noise component. However, this study further assumes that individuals' transformed efficiency also depends on certain characteristics. The extended spatial autoregressive model for the transformed technical efficiency s t is written as:

ARTICLE IN PRESS
with its corresponding Data Generating Process (DGP) being: where I is the N × N identity matrix, Z t is a N × L matrix of L farm characteristics observed in time t that may affect transformed technical efficiency s t , δ is the associated L × 1 vector of parameters to be estimated and ξ t is a N × 1 vector of two-sided error terms that are assumed to be Normally distributed with a mean value of 0 and variance σ 2 ξ I . With respect to the model's spatial component in Eq.
(2) , W is the N × N spatial weights matrix and the product Ws t is the so-called spatial lag. For each individual, this component represents a linear combination of transformed efficiency values of neighbors. To achieve this, elements are placed on the spatial weights matrix W such that, for each individual, the spatial lag is a scalar that represents a linear combination of transformed efficiency values taken by neighbors. Finally, the ρ parameter is the spatial autoregressive parameter that scales the spatial lag.
According to Anselin (1988) , LeSage and Pace (2009) and Elhorst (2014) , statements regarding the strength of spatial effects are not directly based on the estimate of the spatial autoregressive parameter ρ from Eq. (2) but on the calculation of the marginal effects of the variables in Z t on (transformed) technical efficiency. These marginal effects are derived based on the DGP presented in Eq. (3) by calculating the derivative of s t with respect to a particular variable contained in Z t . Given that the term (I − ρW ) −1 is involved in the DGP of Eq. (3) , this marginal effect not only provides information on the change of individuals' transformed efficiency due to a change in their characteristic (due to the identity matrix I ) but also due to a (cumulative) change in a neighbors' characteristic (due to the component ρW ). Given the above-presented way of performing inference in a spatial autoregressive model, it becomes obvious that all the previous spatial efficiency studies of Areal et al. (2012) , Fusco and Vidoli (2013) , Tsionas and Michaelides (2016) , Carvalho (2018) and Pede et al. (2018) are not able to draw such inference because they do not include environmental variables as potential determinants of efficiency 3 . Instead, they report the value of ρ but they can not make any statement regarding spatial spillovers and their magnitude, which can be inferred based on the effect of neighbors' characteristics on individuals' efficiency.

Endogeneity & inference
Although not discussed by previous spatial efficiency studies, the model presented in Eq. (2) suffers from endogeneity bias ( Anselin, 1988;LeSage & Pace, 2009 ). This becomes obvious if one replaces s t on the right-hand side of Eq. (2) with s t = ρWs t + Z t δ + ξ t . This yields: Performing the multiplications leads to: In the last expression, the matrix W 2 reflects second-order neighbors (i.e. those that are neighbors to the first-order neighbors). Since the neighbor of the neighbor (i.e. second-order neighbor) of a particular observation includes the observation itself, W 2 has positive elements on the diagonal. Hence, for a particular observation, s t appears both on the left and on the right-hand-side of the equation and the model suffers from simultaneity (see LeSage & Pace, 2009 , page 14).
To overcome such a bias, the current study specifies a lag spatial lag component: This time, if one replaces s t−1 on the right-hand side of Eq. (6) with s t−1 = ρWs t−2 + Z δ + ξ t−1 yields: Performing the multiplications in the above equation leads to: In the above case, W 2 has again positive elements on the diagonal entries as the second-order neighbors of a particular observation include the observation itself. However, for a particular observation, while s t appears on the left-hand-side, the right-hand-side contains s t−2 and the endogeneity issue is wiped out.
The disadvantage of the above procedure though, is that one is not able to write the DGP of the model presented in Eq. (6) anymore as this is done in Eq. (3) . This is because when the lag spatial lag component is moved to the left-hand side of Eq. (6) , one is not able to take a common factor the dependent variable since s t = s t−1 . Not being able to write the DGP process entails that one is not able to draw the typical inference for a spatial autoregressive model and (qualitative) statements can only be made based on the spatial autoregressive parameter ρ.
Nevertheless, this study argues that given the modelling approach presented in Eq. (6) , the typical inference for the spatial autoregressive model can still be drawn but in a long-run framework. This can be achieved if one applies recursive substitution for past values of s over t periods in Eq. (6) and assumes stationarity 4 . Following these steps, the model that initializes the stochastic process of s and can be viewed as a long-run equilibrium is: where s 0 stands for the N × 1 vector of long-run transformed efficiency and ξ 0 is a N × 1 vector of two-sided long-run error terms assumed to be Normally distributed with mean 0 and variance σ 2 ξ 0 I . Obviously, the model in Eq. (9) is the same as the DGP of a typical spatial autoregressive model (as this is presented in Eq. (3) ) with the difference being that it now concerns long-run transformed efficiency. Given the transformation used for the efficiency component, the long-run technical efficiency ( TE 0 ) corresponds to the expected value of Accordingly, the derivative of the expected long-run technical efficiency with respect to the l th explanatory variable contained in Z is calculated as 5 : [m5G;November 14, 2019;9:28 ] Evaluation of the above presented derivative yields a N × N matrix from which three different types of marginal effects can be derived. First, each diagonal entry represents the so-called direct marginal effect that reveals the percentage change in the longrun efficiency of individual i as a result of a 1 unit change in an explanatory variable. Second, the row sum for each individual i (excluding the respective diagonal element) yields the socalled indirect marginal effect that captures the percentage change in the long-run efficiency of individual i due to a 1 unit (cumulative) change in the explanatory variable of neighbors. Third, adding the direct and indirect marginal effects yields the so-called total marginal effect. Given that each individual possesses a marginal effect, LeSage and Pace (2009) proposes calculation of average marginal effects. Specifically, the mean of the diagonal elements represents the average direct marginal effect, the mean of row sums yields the average indirect marginal effect, while adding these two yields the average total marginal effect.

Bayesian estimation & simulation
The model in Eqs. (1) , (6) and (9) is estimated using Bayesian techniques. For convenience, all parameters to be estimated are In the equations that follow, the latent data a and s t are placed inside curly brackets and are obtained as by-products of the utilized sampler. The complete data likelihood of the model is written as: The first line of Eq. (11) is the probability density function (pdf) of the Normal distribution imposed on the error term v t from Eq. (1) , the second line is the pdf of the Normal distribution imposed on the error term ξ 0 from Eq. (9) , the third line is the pdf of the Normal distribution imposed on the error term ξ t from Eq. (6) and the fourth line is the pdf of the Normal distribution imposed on the error term a from Eq. (1) . Using Bayes' rule, the joint posterior density of the model's parameters contained in θ and the latent variables a and s t is proportional to the product of the complete data likelihood and the priors imposed on the parameters: The combination of data augmentation with Markov Chain Monte Carlo (MCMC) simulation is used to integrate the latent variables a and s t from the likelihood and draw samples from the posterior 6 . The performance of the proposed estimation method is evaluated using simulated data in an attempt to verify that the utilized lag spatial autoregressive model does not yield any biases. For this purpose, the following strategy is adopted: • A panel dataset of 100 individuals with 10 time observations per individual is constructed. 6 A technical appendix that presents the model setup and all related Bayesian quantities including the likelihood, the priors, the posterior and the full conditionals of the parameters is available in the supplementary materials. • Three X t variables and two Z variables are constructed as random draws from standard normal distributions. • Latitude and longitude data are obtained as random draws from a uniform distribution and are used to calculate the distance between individuals. Following Wang, Kockelman, and Wang (2013) , the minimum and maximum values of the unifrom distribution are set equal to those of the real dataset to offer a somewhat realistic geographic setting. The spatial weights matrix W is then constructed based on the inverse distances, with zeros being specified on the diagonal and in the entries where the distance is above the minimum value at which all individuals have at least one neighbor (as this is the strategy followed in the empirical application with more details being provided in the following section). • Data on s t and s 0 are constructed according to Eq. (6) and Eq. (9) , respectively. • Data on the dependent variable y t are generated according to Eq.
The sampling scheme used in this simulation study involves the following: (a) a burn-in phase of 10,0 0 0 iterations is used to remove the influence of the sampler's initial values, (b) 20 0,0 0 0 iterations are run and (c) one out of two iterations are discarded to remove potential autocorrelations 7 . The results from the simulation study are presented in Table 1 . As it is obvious, all parameters are well identified as the mean values are virtually the same as the true ones. Focusing on the parameter of interest -the spatial autoregressive parameter ρ -the estimation yields a mean value of 0.306, which is virtually the same as the true value (i.e. 0.300). This result adds credibility to the assumption that the lag spatial autoregressive model utilized in this study does not suffer from biases.

Dataset & model specification
The utilized dataset stems from the Dutch FADN system as collected by Wageningen Economic Research and contains 1552 observations of 194 specialized Dutch dairy farms that cover the period 2009-2016 8 . The production frontier stacked over individuals, is specified as translog in inputs and time trend for each time JID: EOR [m5G;November 14, 2019;9:28 ] period:

ARTICLE IN PRESS
Output consists of three subcategories: milk and milk products, turnover and growth of cattle, and crop and other products. A Törnqvist index was constructed using price indices from Eurostat for each category. These price indices are 'Milk', 'Animals' and 'Crop output'. The total reported value of output was then deflated using the Törnqvist index. Capital is composed of the subcategories buildings and machinery. A Törnqvist index was again constructed using the Eurostat price indices 'Buildings' and 'Machinery and other equipment'. The total reported value of capital was then deflated using the Törnqvist index. Labor and land represent the total working hours and the total hectares of utilized agricultural area respectively, while animals consist of the total number of livestock units 9 .  (6) and Eq. (9) contains the following variables: the subsidies that farms receive measured in € 10,0 0 0, the economic size of farms measured in 10 0 European Size Units (ESU) and farmers' degree of specialization measured as milk revenues over total revenues 10 . The choice of the aforementioned variables as potential determinants of (transformed) efficiency is based on the following reasoning: subsidies can either be linked to higher efficiency if used for investing in new technologies ( Kumbhakar & Bokusheva, 2009 ) or lower efficiency if viewed as a source of income that lowers the working motivation ( Zhu, Demeter, & Oude Lansink, 2012 ). Economic farm size is expected to be positively associated with efficiency because it captures the impact of economies of scale, which can affect efficiency positively ( Zhu et al., 2012 ). Finally, specialized farmers are expected to be more efficient because of the higher experience that they gain when focusing on one production activity and the associated lower probability of misusing their resources ( Zhu et al., 2012 ). Summary statistics of all the above specified variables appear in Table 2 . Prior to estimation, the output and inputs are normalized by their respective geometric means, while the trend variable is normalized by its arithmetic mean. The normalization used makes 9 The empirical literature mostly treats animals as an input ( Areal et al., 2012;Dakpo, Jeanneaux, & Latruffe, 2017;Emvalomatis et al., 2011;Kumbhakar, Ghosh, & McGuckin, 1991;Njuki, Bravo-Ureta, & Mukherjee, 2016;Skevas et al., 2018b ). This is because animals are used to produce milk, meat and even more animals. The last is captured by including growth of cattle as an output and therefore animals are specified as an input. Furthermore, in the results that follow, the MCMC chain of the parameter 'animals' is fully converged, while formal model comparison (based on Bayes factors) between a model that excludes and a model that includes 'animals' as an input suggests that the latter is favored by the data. Model comparison quantities and the figure of the 'animals' MCMC chain can be provided upon request. 10 Variation of the Z variables over time is negligible. The average coefficient of variation is 0.102 for subsidies, 0.152 for ESU and 0.034 for specialization. This adds further credibility to the specification of time-invariant Z variables in Eqs. (6) and (9) . the first-order term parameters interpretable as output elasticities evaluated at the geometric mean of the data.

Specification of the spatial weights matrix
The spatial weights matrix W is constructed based on the distance between farms, which is calculated using the information that the utilized dataset provides on farms' latitude and longitude. Recognizing that individuals may be influenced more by closer than distant neighbors, an inverse distance matrix is specified with entries w i j = 1 d i j , where d ij is the distance between individual i and j . The typical approach of setting a distance cut-off point d * is followed and all spatial weights w i j outside this distance are zero. Instead of setting an arbitrary threshold, this study chooses the threshold of 18km, which corresponds to the minimum distance at which all farms in the sample have at least one neighbor. This is a standard approach followed by several studies in the spatial econometrics literature including Osland (2010) and Marasteanu and Jaenicke (2016) . Furthermore, all diagonal elements w ii are also set equal to zero so that individuals are not defined as neighbors to themselves. Finally, the conventional practice of standardizing the spatial weights matrix W is followed by scaling all elements by its maximum eigenvalue as in Vega and Elhorst (2015) .

Specification of priors
The following priors are imposed on the parameters to be estimated: • A multivariate normal density is used for the joint prior densities of β and δ. In both cases, the prior means are set equal to zero and the prior covariance matrices take the value of 10 0 0 on their diagonal entries. Setting such high prior variances implies that the prior beliefs will have a negligible impact on the results.

Results
Following the same sampling scheme used in the simulation study yields the posterior moments of the production frontier parameters that are presented in Table 3 11 , 12 . As mentioned before, the geometric mean transformation used for outputs and inputs makes the first-order parameters with respect to inputs interpretable as output elasticities evaluated at the geometric mean of the data. All output elasticities are positive, which implies that the monotonicity condition (i.e. non-negative elasticities) of the production function is fulfilled at the geometric mean of the data. Apart from the elasticity with respect to capital, the rest of the elasticities are statistically significant since their corresponding 95% credible intervals do not contain zero. The reported output 11 As mentioned in the previous section, all reported results are based on the distance cut-off point of 18km, which corresponds to the minimum distance at which all farms in the sample have at least one neighbor. Robustness checks with respect to higher arbitrary cut-off points of 20km and 22km are performed yielding virtually the same coefficient estimates. The results of these robustness checks can be provided upon request. 12 A production frontier is estimated because although specialized Dutch dairy farms produce multiple outputs, the share of outputs other than milk to the total output is very low as this is manifested by the mean value of specialization in milk production presented in Table 2 . Besides, robustness checks reveal that separating other outputs from milk and employing a distance function does not result to changes in the estimated elasticities and the efficiency scores. Furthermore, the production frontier specification is favored by the data when compared to a distance function specification according to Bayes factors. The results of the robustness checks can be provided upon request. elasticities have also magnitudes that are consistent with the findings of previous studies that examined Dutch and other western European countries' dairy farms efficiency. Specifically, animals, intermediate inputs and feed have the highest effect on production, which is also reported in the studies of Emvalomatis et al. (2011) andSkevas, Zhu, Shestalova, andEmvalomatis (2018d) for the case of Dutch dairy farms. Additionally, Areal et al. (2012) find that UK's dairy farms output is mainly affected by animals, while Skevas, Emvalomatis, and Brümmer (2018a) and Skevas, Emvalomatis, and Brümmer (2018c) also find a strong association between output and animals, feed and intermediate inputs for the case of German dairy farms. Adding the output elasticities yields a scale elasticity of 0.9. Dutch dairy farms operate under the decreasing returns to scale part of the technology with a probability of 99.6%. Furthermore, farms do not experience a significant technological regress as the estimated parameter with respect to trend is negative but insignificant.
Average technical efficiency (obtained as the average of exp (s t ) / (1 + exp (s t )) across all farms and years is 0.843. That is, Dutch dairy farms produce 84.3% of what is feasible given the amount of inputs used. Average long-run technical efficiency across all farms is 0.845. The fact that average efficiency is almost identical with average long-run efficiency implies that the time-span captured by the utilized data is close to its equilibrium. The reported average efficiency is similar to the mean value of 83.1% reported by Emvalomatis et al. (2011) for the case of Dutch dairy farms and to the average score of 84% found by Areal et al. (2012) for UK dairy farms. However, Skevas et al. (2018c) reported a mean efficiency estimate of 70% for German dairy farms, while Dakpo et al. (2017) reported an average efficiency score of 74% for the case of French dairy farms. These discrepancies may arise because Germany and France are larger countries with more heterogeneous farms that can lead to lower average efficiency estimates. Additionally, the discrepancies in the mean efficiency scores may be attributed to modelling differences as Skevas et al. (2018c) estimated a non-spatial dynamic hierarchical model, while Dakpo et al. (2017) also accounted for undesirable outputs.
Moving to the estimates with respect to the determinants of transformed technical efficiency s t , Table 4 presents their posterior moments. The estimate of the spatial autoregressive parameter ρ points towards the existence of spatial dependence in farmers' transformed efficiency levels. Subsidies and specialization have a positive and significant impact on transformed efficiency, while the coefficient with respect to ESU is statistically insignificant.
As consistently mentioned in this article, although the estimate of the spatial autoregressive parameter ρ provides an indication of spatial dependence in farmers' transformed efficiency scores, one is unable to make statements regarding the nature and magnitude of spatial effects. Such statements are only possible through the derivation of the (average) direct and indirect marginal effects of the utilized farm characteristics on transformed efficiency, with the latter providing evidence for the type of spatial spillovers and their magnitude. Given that the modelling approach used makes it possible to draw inference in a long-run framework and that interest lies on efficiency rather on its transformed quantity, the (average) direct and indirect marginal effects of the utilized JID: EOR [m5G;November 14, 2019;9:28 ]  farm characteristics on long-run technical efficiency are derived according to Eq. (10) and are presented in Table 5 13 14 .

ARTICLE IN PRESS
The average direct marginal effect of subsidies is 0.037. That is, a 1 unit (i.e. € 10,0 0 0) increase in farmers' subsidies is associated with 3.7% higher long-run technical efficiency, on average. A possible explanation for this finding is that subsidies act as a source of credit that is used for investment purposes, thus increasing farms' long-run efficiency ( Kumbhakar & Bokusheva, 2009 ). The same finding is reported by Rizov, Pokrivcak, and Ciaian (2013) for Danish dairy farms, while Zhu et al. (2012) report a negative association between subsidies and efficiency for the case of Dutch dairy farms. An explanation for this contradicting finding is that the time period studied by Zhu et al. (2012) covered the beginnings of the decoupling reform of the Common Agricultural Policy (CAP) in 2003, which offered a strong incentive to Dutch farmers to rely on subsidies and refrain from productive activities. However, the more recent time period considered in the current study covers the new 2013 CAP reform in which subsidies are also provided to encourage farmers to use more up-to-date and environmentally friendly technologies that can increase their efficiency. The average indirect marginal effect of subsidies is 0.019. This means that, on average, a 1 unit cumulative increase in neighbors' subsidies is linked to a 1.9% higher long-run technical efficiency for individuals. Therefore, a positive spillover effect arises when neighbors' budget increases and allows them to invest in up-to-date technologies. A possible explanation for this finding is that farmers who are surrounded by neighbors that excel using more efficient production technologies, become more motivated to perform better in order to remain competitive and avoid a potential market exit.
The average direct and indirect marginal effects of ESU do not significantly affect long-run technical efficiency. In terms of specialization, the corresponding average direct marginal effect is estimated at 0.022. This means that on average, a 1% increase in individuals' degree of specialization is associated with 2.2% higher long-run technical efficiency. This result may be explained by the fact that focusing on a particular production activity makes farmers more experienced and less prone in committing production mistakes ( Zhu et al., 2012 ). The same finding was reported by Latruffe, Davidova, and Balcombe (2008) and Sauer and Latacz-Lohmann (2015) for the case of Polish and German dairy farms, respectively. The average indirect marginal effect of specialization is 0.011. This suggests that a 1% cumulative increase in neighbors' levels of specialization is linked to a 1.1% higher long-run technical efficiency for individuals, on average. Such a result highlights a potential 13 Note that normally, the estimated coefficients of the variables included in the spatial autoregressive model that are presented in Table 4 should be close in magnitude to the derived direct marginal effects presented in Table 5 . However, the coefficient estimates deviate from each other because Table 4 reports the effects of the utilized variables on transformed technical efficiency s t , while Table 5 presents the corresponding direct marginal effects on long-run efficiency TE 0 . Instead, the magnitudes of the coefficient estimates of the variables in the spatial autoregressive model would be similar to the respective direct marginal effects if the latter would be calculated for long-run transformed efficiency s 0 , which is not however a standard practice in the literature as interest mostly lies on making quantitative statements for efficiency. 14 Table 5 presents the mean values of the marginal effects and omits their standard deviation and 95% credible intervals for space-saving purposes. Note, however, that significance of the estimates is the same as in the case of transformed efficiency. learning spillover effect through which farmers learn how to use their resources efficiently from their more experienced neighbors, which results in their efficiency gains in the long-run.

Conclusions
This article proposes an extension to the typical spatial autoregressive efficiency model that concerns the inclusion of firm characteristics as further potential determinants of efficiency. Such an extension provides theoretical as well as methodological advantages. Theoretically-wise, the restrictive assumption of previous spatial efficiency studies that individuals' efficiency is only related to their neighbors' efficiency and a noise component is relaxed by further assuming that it also depends on certain (widely studied) firm characteristics. Methodologically-wise, past spatial efficiency studies are unable to explore the nature and magnitude of spatial spillovers because they inevitably base their inference only on the estimated parameter associated with neighbors' efficiency, which merely allows to draw qualitative conclusions with regard to spatial effects ( Anselin, 1988;Elhorst, 2014;LeSage & Pace, 2009 ). In contrast, accounting for firm characteristics allows this study to perform the standard inference presented in spatial autoregressive models that concerns the derivation of direct and indirect marginal effects. The former reveal changes in individuals' efficiency as a result of changes in their own characteristics, while the latter quantify changes in individuals' efficiency due to changes in their neighbors' characteristics. Furthermore, unlike previous studies, the present one tackles the issue of endogeneity in the utilized spatial autoregressive efficiency model, which arises due to simultaneity, as the dependent variable for any individual appears both on the left and on the right-hand-side of the model. A lag spatial lag component is specified that wipes out endogeneity and makes inference to be performed in a long-run equilibrium framework. Finally, all the above are supported by the use of a very recent panel dataset for specialized Dutch dairy farms that contains exact information of their location based on latitude and longitude.
A production frontier is estimated while the extended spatial autoregressive model is imposed on transformed efficiency. Average technical efficiency across all individuals and years is estimated at 84.3% and average long-run technical efficiency at 84.5%. This result suggests that Dutch dairy farms are close to their equilibrium while highlighting their potential to further improve their efficiency levels. The efficiency scores exhibit spatial dependence, which is however a finding that does not directly provide any information on the type and magnitude of spatial effects. Inf erence in the estimated spatial autoregressive model is rather based on the derivation of direct and indirect marginal effects of the utilized farm characteristics on long-run technical efficiency, with the latter providing evidence on the nature and strength of spatial spillovers. With respect to subsidies, an increase in farmers' subsidies is associated with higher long-run efficiency revealing that subsidies can ameliorate farms' financial constraints enabling them to adopt new technologies that increase their long-run efficiency. An increase in neighbors' subsidies is also linked to higher long-run efficiency for adjoining farmers. This finding suggests that when neighboring farmers' financial capacity is strengthened through the provision of subsidies that allows them to invest in new technologies and increase their long-run efficiency, this can motivate adjoining farmers to increase their long-run efficiency in order to keep up with their neighbors and remain competitive. In terms of specialization, when farmers become more specialized in milk production, higher long-run efficiency levels are attained, which can be due to the higher experience that makes them less prone in committing production mistakes. An increase in neighbors' levels of specialization, is associated with higher long-run efficiency for adjoining farmers. This result highlights a potential learning spillover effect, which can arise through communication/imitation between neighboring farmers. That is, farmers increase their long-run efficiency by learning how to use their resources more efficiently by more experienced neighbors.
Finally, several policy conclusions can be drawn from the results of the present study. First of all, the fact that feed and the volume of animals play the most important role in Dutch dairy farms' production provides a signal that any changes in the use of these inputs can significantly alter production. Therefore, policymakers need to protect farmers if experiencing high feed prices as this may force them to lower their feed use and consequently their output. Furthermore, policies aiming at the extensification of the sector by reducing herd sizes, should take into account this large association between the number of animals and production in order to prevent significant output losses. Additionally, and as mentioned above, there is a scope for Dutch dairy farms to further improve their efficiency levels. According to the empirical findings of this study, this can be materialized by granting subsidies to farmers (provided that they are used for investment purposes) and by promoting specialization in single production activities. Furthermore, efficiency increases can be realized by taking into account the existing motivation and learning spillovers between neighboring farmers that are reported in the present study. For instance, policy-makers can take advantage of the communication between neighboring farmers and promote more efficiently the adoption of state-of-the-art technologies that can increase farms' efficiency. This can be achieved by strengthening local farmers associations and intensifying farm visits, which can facilitate the dissemination of information regarding the existence of new technologies.

Disclaimer
The Data used in the present work stem from the Dutch FADN system as collected by Wageningen Economic Research. The Centre of Economic Information (CEI) has provided access to these data. Results shown are and remain entirely the responsibility of the author; neither they represent Wageningen Economic Research/CEI views nor constitute official statistics.

Appendix
This Appendix presents the derivation of the expected long-run technical efficiency ( TE 0 ) as well as the calculation of the corresponding marginal effects (i.e. the derivative of TE 0 with respect to the l th variable contained in Z ).
To derive the expected value of TE 0 one needs to start from the spatial autoregressive equation: A recursive substitution for the above equation gives: Then, the derivative of the expected long-run technical efficiency TE 0 with respect to the l th variable contained in Z can be calculated using the quotient rule: which after further manipulations results in:

Supplementary material
Supplementary material associated with this article can be found, in the online version, at doi: 10.1016/j.ejor.2019.10.033 .