ENDOGENEITY IN STOCHASTIC PRODUCTION FRONTIER WITH ONE AND TWO-STEP MODELS: AN APPLICATION WITH MUNICIPAL DATA FROM THE BRAZILIAN AGRICULTURAL CENSUS

ABSTRACT

Stochastic production frontier models are widely used in microeconometrics and, in the last decades, have been proven to be versatile in their range of applications. However, there are few studies concerning endogeneity in stochastic production frontier models. Here we present two stochastic production frontier models with endogenous variables based on the main distributions for the technical inefficiency. We also derive analytic gradient vectors to obtain the best performance at a reasonable computational time cost. The methodology presented here is based on one and two-step maximum likelihood estimation, allows for endogeneity and heteroscedasticity in relation to one or both error terms, and is implemented in R language. Finally, we illustrate an application with municipal data from the Brazilian agricultural census. The results show that capital dominates the production function, credit access and technical assistance are endogenous, and income concentration seems to impede productive inclusion through the more intensive use of technology.

Keywords:
endogeneity; maximum likelihood method; stochastic production frontier

1 INTRODUCTION

Stochastic production frontier models were simultaneously introduced by Aigner et al. (1977); Meeusen & van Den Broeck (1977) and are widely used in microeconometrics. A stochastic frontier model is a random-effects model designed to estimate the technical efficiency of decision-making units or producers through production or cost functions. The literature includes many empirical examples from various fields, such as agriculture, banking, and health.

By using cross-sectional data on the quantities of k inputs to produce a single product for each of n producers, we can write a stochastic production frontier model as

$y_i=f(x_i;\beta )\exp \left(v_i\right){\mathcal{E}}_i=f(x_i;\beta )\exp (v_i-u_i),e_i=v_i-u_i,i=1,2,\dots ,n,$ (1)

where y _i is the dependent variable, x _i is a vector of k inputs used by the i-th producer, β is a $k\times 1$ vector of technology parameters to be estimated, and $f(x_i;\beta )\exp \left(v_i\right)$ is the stochastic production frontier, which consists of two parts: a deterministic part, $f(x_i;\beta )$, common to all producers, and a producer-specific part, exp(v _i ), which captures the effect of each producer’s specific random shocks.

Since the component ${\mathcal{E}}_i=\exp (-u_i)$ is the production-oriented technical efficiency of the i-th producer, we have

${\mathcal{E}}_i=\frac{y_i}{f(x_i;\beta )\exp \left(v_i\right)},0<{\mathcal{E}}_i<1,$ (2)

which defines technical efficiency as the ratio between observed production and maximum feasible production, that is, ℰ _i provides a measure of the observed production deficit of each producer relative to the maximum feasible production in an environment characterized by exp(v _i ) (Kumbhakar & Lovell, 2003).

Estimates of each producer’s technical efficiency depend on decomposing e _i and are typically derived from the conditional expectation of exp(−u _i ) given e _i , which vary according to the probability density functions of both v _i and u _i .

Consider the stochastic production frontier model of (1) in the log-linear Cobb-Douglas form,

$\ln y_i={\beta }_0+\sum _{j=1}^k{\beta }_j\ln x_{ji}+v_i-u_i.$ (3)

In this specification, $\ln f(x_i;\beta ).$ produces a linear model on β, with the usual Gaussian noise, $v~\mathrm{\mathbb{N}}(0,{\sigma }_v^2)$, and a random effect, $u~f_u$, representing the unit’s technical inefficiency.

In empirical studies, it is common to assume that u _i has a half-normal istribution. However, other distribution assumptions concerning the one-sided error term (u _i ) have been proposed, such as the exponential and the truncated normal distribution. We can use maximum likelihood methods to estimate the parameters of such models. In this context, the inefficiency term is a latent variable that must be integrated when calculating the likelihood. Depending on the choice of the u _i density function, the likelihood calculation will require numerical integration. Greene (1990) and Andrade & Souza (2017) discuss approximation techniques and their accuracy. Consequently, the normal/gamma model is not available in most statistical and econometric tools for stochastic frontier analysis. Alternatively, motivated by the hierarchical structure of the stochastic frontier models, Andrade & Souza (2019) use the Expectation-Maximization (EM) algorithm to perform stochastic frontier analysis. In their work, the Expectation-Maximization calculations resulted in simple algorithms with closed-form expressions for the half-normal and exponential models and more elaborate versions for the truncated normal and gamma models.

Moreover, it is common to assume a dependency of the inefficiency component on several factors that may influence performance. With the half-normal and the exponential distributions, this is achieved by postulating a dependency of the log-likelihood on a linear construct (ζ 'z), defined by a set of contextual variables z. For the truncated normal, the mean of the underlying normal is defined by the linear construct. Often, we introduce heteroscedasticity in the exponential and half-normal cases, assuming that the Gaussian error terms may not have constant variance. This is achieved by also assuming that ${\sigma }_v^2$ is dependent on a linear function of a set of, known covariates w, i.e., ${\sigma }_v^2=\exp \{\rho \text{'}w\}$.

Another critical issue, as pointed out by Cazals et al. (2016), is that there may be factors (observable for the firm but unobservable for the econometrician) that affect both the choice of regressors and inefficiency levels. These factors would be endogenous variables. Consequently, failure to recognize endogeneity in the production function and the inefficiency components will lead to inconsistencies in the estimation process.

Production frontier analysis aims to identify best production practices and the importance of external factors, endogenous or not, affecting the production function and the technical efficiency component. In particular, we are interested in identifying the effect on the production of variables related to market imperfections. Following Souza et al. (2017) and Souza & Gomes (2018), market imperfections occur when farmers are subjected to different market conditions depending on their income. Large-scale farmers generally access lower input prices and sell their products at lower prices, making competition harder for small farmers. Therefore, market imperfections are typically associated with infrastructure, environment control requirements, and the presence of technical assistance. Identifying these factors and estimating the corresponding elasticities are fundamental for public policies envisaging productive inclusion.

Some programs have stochastic frontier analysis implementations available. However, often the error terms of stochastic frontier models do not have constant variance. Most stochastic frontier packages or routines available do not allow fit models considering heteroscedastic error components. For example, the sfa package in the R language permits specify half-normal, exponential, or truncated normal distributions for the one-sided error term, u. However, it does not allow you to insert a set of covariates to model the error terms. Besides that, many of them only deal with the exogeneity assumption, failing to address the endogeneity that may exist when one or more frontier or inefficiency variables correlate with the two-sided error term, v _i. Such limitations make it of interest to provide routines that enable modeling under these scenarios.

Additionally, many programs use numerical procedures to estimate parameters. However, by including the analytic gradient vector in the estimation process, we can considerably improve the convergence rate of the iterative optimization procedure. Unfortunately, in the stochastic frontier analysis literature, no work is yet available to illustrate the expressions of these gradients when endogenous variables are present.

Therefore, this paper aims to implement the prediction of the producers’ technical efficiency level from stochastic production frontier models containing endogenous and exogenous variables through one and two-step maximum likelihood estimation procedure, based on Karakaplan & Kutlu (2015), considering the main specifications for the inefficiency term (half-norm, exponential, and truncated normal distributions), as well as to derive the analytic gradients of these models. Furthermore, we implemented a routine in the R language for these models, which makes it possible to deal with endogeneity and heteroscedasticity in any term of the model.

This article begins with a brief review of stochastic frontier models. Then, in Section 2 we cover the stochastic production frontier literature in the presence of endogeneity; Section 3 describes the database used in the application; Section 4 concerns the models and gradients derived; Section 5 discusses the results of applying these techniques to real data. Finally, Section 6 presents the conclusions reached and possible further studies.

2 LITERATURE

As the consistency of the usual stochastic frontier estimators depends on the exogeneity of the regressors, standard estimators do not deal with the endogeneity that exists if the frontier or inefficiency variables correlate with the two-sided error term, v _i , leading to inconsistent parameter estimation. Mutter et al. (2013) explains why omitting the variable causing endogeneity is not a feasible solution. Consequently, dealing with endogeneity in the stochastic frontier analysis is relatively more complicated than in standard regression models due to the unique nature of their error terms.

Some of the first stochastic frontier articles to tackle endogeneity are Guan et al. (2009); Kutlu (2010); Tran & Tsionas (2013). In these studies, variables can be endogenous because they correlate with v _i but not with u _i . Guan et al. (2009) propose a two-step estimation method to handle the endogenous regressors in the model. In the first step, they obtain consistent estimates of the frontier parameters by a generalized method of moments. In the second step, they use the residuals from the first step as a dependent variable and then use the maximum likelihood method to estimate excess capital capacity. Kutlu (2010) describes a model for dealing with endogeneity by maximum likelihood estimation at one and two-step, where he estimates the time-varying technical efficiency in the presence of endogenous regressors by using a modified version of the Battese & Coelli (1992) estimator. Tran & Tsionas (2013) propose a variation of Kutlu (2010) by a generalized method of moments. However, these model assumptions are insufficient to tackle endogeneity due to the correlation between the error terms, u _i , and v _i .

Tran & Tsionas (2015) and Amsler et al. (2016) handle endogeneity with a copula approach. Tran & Tsionas (2015) use a copula function to directly model the correlation between the endogenous regressors and the composed error term. While Amsler et al. (2016) allows the endogeneity of the regressors concerning statistical noise and inefficiency separately. A copula approach allows more general correlation structures when modeling endogeneity. However, this method is computationally intensive and requires choosing a suitable copula. In addition, the models proposed by Tran & Tsionas (2015); Amsler et al. (2016), as well as Guan et al. (2009); Kutlu (2010); Tran & Tsionas (2013) do not allow contextual variables that affect inefficiency.

Griffiths & Hajargasht (2016); Karakaplan & Kutlu (2015) handle endogeneity concerning one and two-sided errors and the correlation between them. In addition to considering environmental variables that affect inefficiency, Griffiths & Hajargasht (2016) presents a Bayesian stochastic frontier model where u _i or v _i or both correlate with the regressors. However, their model is very different from the model proposed by Karakaplan & Kutlu (2015). Instead, Karakaplan & Kutlu (2015) suggest using instrumental variables and a methodology based on a one-step maximum likelihood estimation to obtain a consistent estimator in the presence of endogeneity due to the correlation between the error terms, allowing v _i and u _i to depend on covariates that shape both distributions. In general, one of the main strengths of this model is that it is easier to apply than the copulas approach or Bayesian models, and it is a direct generalization of one of the most used stochastic frontier models - estimators of the Battese & Coelli (1995) type. In their model, Karakaplan & Kutlu (2015) assumes a linear regression by instrumental variables, but the idea can easily be generalized for a nonlinear specification. In addition, they consider a half- normal distribution for u _i . Karakaplan (2017) provides the sfkk module on Stata for this model specification.

Other recent works dealing with endogeneity are Prokhorov et al. (2020); Tsionas et al. (2021); Kumbhakar et al. (2020). Prokhorov et al. (2020) consider the problem of estimating a non- parametric stochastic frontier model with shape restrictions and when some or all regressors are endogenous. They discuss three estimation approaches based on constructing a likelihood with unknown components. Tsionas et al. (2021), using US banking data, propose a Bayesian approach for inference in the stochastic ray production frontier, which can model multiple-input - multiple-output production technologies even in case of zero output quantities. Finally, Kumbhakar et al. (2020) discuss the range of methods developed over the last four decades concerning stochastic frontier analysis.

3 DATA

For the application, we used cross-sectional data from the 2006 Brazilian agricultural census aggregated at the municipal level; from the 2010 Brazilian demographic census; from the National Institute of Research and Educational Studies (INEP), referring to education in 2009; and from the Ministry of Health in 2011. These data are valid for 4965 municipalities, which account for almost 90% of the total number of Brazilian municipalities.

The production model assumes as a dependent variable the gross income of the rural establishments in reais (income), i.e., the total value of the agricultural production of the establishments and, as inputs, the expenses on land (land), labor (labor) and capital (techinputs, technological inputs). These variables were extracted from the 2006 Brazilian agricultural census database and aggregated at the municipal level.

The contextual variables affecting production - keys to balance market imperfections - are aggregate indicators referring to the social (social), demographic (demographic) and environmental (environment) characteristics of the rural development. We also considered variables related to credit access (financing, total financing per farm), technical assistance (techassist, proportion of farms who received technical assistance), an indicator of income concentration per municipality (gini) and dummies for regional effects (regions, variables indicating county regions).

Except for the region, the other contextual variables were ranked and normalized by the maximum value. This approach lends nonparametric statistical properties to the analysis and circumvents outliers and heteroscedasticity problems.

For the production function, the logarithm of the gross income variable is considered a response variable. As explanatory variables we have the logarithm of the production factors - land, labor and techinputs - and regional dummies. The techassist, financing and gini variables are used to model the ${\sigma }_{ui}^2$ function of the producers’ technical inefficiency level. The social, demographic and environment variables are external instrumental variables in this analysis. We assumed that the techassist and financing variables are potentially endogenous. Both are complex variables that can involve many factors related to the structure of the production unit and are strong candidates for endogeneity.

4 METHODOLOGY

Consider the following stochastic production frontier model with endogenous variables proposed by Karakaplan & Kutlu (2015):

$$\begin{gathered} y_i=x_{1i}^{\top }\beta +v_i-u_i, \\ {x}_i=Z_i\delta +{\epsilon }_i, \\ \left[\begin{array}{c}{\tilde{\epsilon }}_i\\ v_i\end{array}\right]\equiv \left[\begin{array}{c}{\Omega }^{-\frac{1}{2}}{\epsilon }_i\\ v_i\end{array}\right]~N\left(\left[\begin{array}{c}0\\ 0\end{array}\right],\left[\begin{array}{cc}{I}_p& {\sigma }_{vi}\rho \\ {\sigma }_{vi}{\rho }^{\top }& {\sigma }_{vi}^2\end{array}\right]\right). \end{gathered}$$ (4)

In these model, y _i is the natural logarithm of the i-th producer; x _1i is a vector of exogenous and endogenous variables; x _i is a $p\times 1$ vector of all endogenous variables (excluding y _i ); $Z_i=I_p\otimes z_i^{\top }$, where z_i is a $q\times 1$ vector of all exogenous and instruments variables; v _i and ε _i are two-sided error terms; $u_i\ge 0$ is a one-sided error term capturing inefficiency; Ω is the variance-covariance matrix of ε _i ; ${\sigma }_{vi}^2$ is the variance of v _i ; ρ is the vector which represents the correlation between ${\tilde{\epsilon }}_i$ and v _i . In this structure, a variable is endogenous if it is not independent of the two-sided error term, v _i .

This model specifications provide a methodology for dealing with endogeneity in stochastic frontier models in a more general setting. The model considers heteroscedasticity in either component of the composed error term, allowing u _i and v _i to be dependent through covariates that shape both distributions.

Based on a Cholesky decomposition of the variance-covariance matrix of $({\tilde{\epsilon }}_i^{\top },v_i{)}^{\top }$, we have:

$\left[\begin{array}{c}{\tilde{\epsilon }}_i\\ v_i\end{array}\right]=\left[\begin{array}{cc}{I}_p& 0\\ {\sigma }_{vi}{\rho }^{\top }& {\sigma }_{vi}\sqrt{1-{\rho }^{\top }\rho }\end{array}\right]\left[\begin{array}{c}{\tilde{\epsilon }}_i\\ {\tilde{w}}_i\end{array}\right],$ (5)

where ${\tilde{\epsilon }}_i$ and ${\tilde{w}}_i~\mathrm{\mathbb{N}}(0,1)$ are independent. Thus, the frontier equation is expressed by:

$$\begin{gathered} y_i=x_{1i}^{\top }\beta +{\sigma }_{vi}{\rho }^{\top }{\tilde{\epsilon }}_i+w_i-u_i, \\ =x_{1i}^{\top }\beta +\frac{{\sigma }_{wi}}{{\sigma }_{cw}}{\eta }^{\top }(x_i-Z_i\delta )+e_i, \end{gathered}$$ (6)

where $e_i=w_i-u_i;w_i={\sigma }_{vi}\sqrt{1-{\rho }^{\top }\rho }{\tilde{w}}_i={\sigma }_{wi}{\tilde{w}}_i;{\sigma }_{cw}>0$ is a function of the constant term of ${\sigma }_{wi};\eta ={\sigma }_{cw}{\Omega }^{-\frac{1}{2}}\rho /\sqrt{1-{\rho }^{\top }\rho }$. Thus, when there is no heteroscedasticity in w _i , σ _wi = σ _cw , so that:

$y_i=x_{1i}^{\top }\beta +{\eta }^{\top }(x_i-Z_i\delta )+e_i.$ (7)

The term e _i is conditionally independent from the regressors given x _i and z _i , and it is possible to directly assume that the conditional distribution of v _i given x _i (and exogenous variables) is a normal distribution with mean equal to $({\sigma }_{wi}/{\sigma }_{cw}){\eta }^{\top }(x_i-Z_i\delta )$. This approach is commonly used to solve the problem of building a consistent estimator in the presence of endogeneity for models with intrinsic nonlinearity such as this model, where $({\sigma }_{wi}/{\sigma }_{cw}){\eta }^{\top }(x_i-Z_i\delta )$ is a bias correction term. Therefore, this approach treats endogeneity as an omitted variable problem.

Let x _2i be a vector of exogenous and endogenous variables and x _3i be a vector of exogenous and endogenous variables, which can share the same variables with x _1i and x _2i . We assume that ${\sigma }_{ui}^2=\exp \left(x_{2i}^{\top }{\phi }_u\right),{\sigma }_{wi}^2=\exp \left(x_{3i}^{\top }{\phi }_w\right)$ and ${\sigma }_{cw}^2=\exp \left({\phi }_{cw}\right)$, where $\phi =({\phi }_u^{\top },{\phi }_w^{\top }{)}^{\top }$ is the vector of parameters capturing heteroscedasticity, and φ _cw is the coefficient of the constant term for $x_{3i}^{\top }{\phi }_w$.

The proposed specifications for the inefficiency term are: half-normal, $u_i~{\mathrm{\mathbb{N}}}^{+}(0,{\sigma }_{ui}^2)$, exponential, $u_i~\mathrm{Exp}\left({\sigma }_{ui}\right)$, or truncated normal, $u_i~{\mathrm{\mathbb{N}}}^{+}({\mu }_i,{\sigma }_u^2)$, with ${\mu }_i=x_{2i}^{\top }\tau$ being the mean of the truncated normal distribution, where τ is the vector of parameters capturing heteroscedasticity. For the truncated normal distribution, it is assumed that only µ _i is a function of covariates, while ${\sigma }_{ui}^2$ and ${\sigma }_{wi}^2$ are constant terms.

The log-likelihood of the stochastic production frontier model decomposes into two parts: $\ln L\left(\theta \right)=\ln L_{y|x}\left(\theta \right)+\ln L_x\left(\theta \right)$,

$$\begin{gathered} \ln L_{y|x}\left(\theta \right)\propto \sum _{i=1}^n\left\{-\frac{1}{2}\ln {\sigma }_i^2+\ln \Phi \left(-\frac{e_i{\lambda }_i}{{\sigma }_i}\right)-\frac{e_i^2}{2{\sigma }_i^2}\right\}, \\ \ln L_{y|x}\left(\theta \right)\propto \sum _{i=1}^n\left\{-\frac{1}{2}\ln {\sigma }_{ui}^2+\frac{{\sigma }_{wi}^2}{2{\sigma }_{ui}^2}+\ln \Phi \left(\frac{-e_i-\frac{{\sigma }_{wi}^2}{{\sigma }_{ui}}}{{\sigma }_{wi}}\right)+\frac{e_i}{{\sigma }_{ui}}\right\}, \\ \ln L_{y|x}\left(\theta \right)\propto \sum _{i=1}^n\left\{-\frac{1}{2}\ln {\sigma }^2-\ln \Phi \left(\frac{{\mu }_i}{{\sigma }_u}\right)+\ln \Phi \left(\frac{{\mu }_i}{\sigma \lambda }-\frac{e_i\lambda }{\sigma }\right)-\frac{1}{2}{\left(\frac{e_i+{\mu }_i}{\sigma }\right)}^2\right\}, \\ \ln L_x\left(\theta \right)=\sum _{i=1}^n\left\{\frac{-p\ln \left(2\pi \right)-\ln \left(\right|\Omega \left|\right)-{\epsilon }_i^{\top }{\Omega }^{-1}{\epsilon }_i}{2}\right\}, \\ {e}_i=y_i-x_{1i}^{\top }\beta -\frac{{\sigma }_{wi}}{{\sigma }_{cw}}{\eta }^{\top }(x_i-Z_i\delta ), \\ {\epsilon }_i=x_i-Z_i\delta , \\ {\sigma }_i^2={\sigma }_{ui}^2+{\sigma }_{wi}^2, \\ {\lambda }_i=\frac{{\sigma }_{ui}}{{\sigma }_{wi}}, \\ {\mu }_i\ne 0\iff u_i~{\mathrm{\mathbb{N}}}^{+}({\mu }_i,{\sigma }_u^2), \end{gathered}$$ (8)

where $y=(y_1,\dots ,y_n{)}^{\top }$ is the vector of dependent variable, $x=(x_1^{\top },\dots ,x_n^{\top }{)}^{\top }$ is the matrix of endogenous variables in the model and $\theta =({\beta }^{\top },{\eta }^{\top },{\phi }^{\top },{\delta }^{\top },{\tau }^{\top }{)}^{\top }$ is the vector of coefficients. Note that x follows a multivariate normal distribution if the number of endogenous variables is greater than one, and univariate normal otherwise. Whereas y|x follows a normal/half-normal, normal/exponential or normal/truncated normal distribution, respectively. ln L _x (θ) is added to ln L _y|x (θ) and e _i is adjusted by the $({\sigma }_{wi}/{\sigma }_{cw}){\eta }^{\top }(x_i-Z_i\delta )$ factor, which solves the problem of inconsistent parameter estimates due to endogenous regressors in x ₁ i and due to the endogenous variables in x _2i . Moreover, it is still possible to test for the presence of endogeneity by testing the null hypothesis that η = 0. More information in Karakaplan (2017).

In this methodology, the parameter vector θ is estimated based on a one-step maximum likelihood estimation method (simultaneously), characterizing a full information maximum likelihood (FIML) approach.

By assuming exponential link function for ${\sigma }_{ui}^2$ and ${\sigma }_{wi}^2$, and ${\mu }_i\in \mathrm{ℝ}$, we obtain the following gradients:

From (8), when assuming a half-normal distribution for u _i , the gradient is

$$\begin{gathered} U\left(\delta \right)=\sum _{i=1}^n{Z}_i^{\top }{\epsilon }_i{\Omega }^{-1}-\sum _{i=1}^n{Z}_i^{\top }\left\{\frac{{\sigma }_{wi}}{{\sigma }_{cw}}\left[\frac{e_i}{{\sigma }_i^2}+\frac{{\lambda }_i}{{\sigma }_i}{A}_i\right]\right\}{\eta }^{\top }, \\ U\left(\beta \right)=\sum _{i=1}^n{x}_{1i}^{\top }\left\{\frac{e_i}{{\sigma }_i^2}+\frac{{\lambda }_i}{{\sigma }_i}{A}_i\right\}, \\ U\left(\eta \right)=\sum _{i=1}^n{\epsilon }_i^{\top }\left\{\frac{{\sigma }_{wi}}{{\sigma }_{cw}}\left[\frac{e_i}{{\sigma }_i^2}+\frac{{\lambda }_i}{{\sigma }_i}{A}_i\right]\right\}, \\ U\left({\phi }_u\right)=\sum _{i=1}^n{x}_{2i}^{\top }\left\{\frac{1}{2{\sigma }_i^2}\left[\frac{e_i^2}{{\sigma }_i^2}-\frac{e_i}{{\lambda }_i{\sigma }_i}{A}_i-1\right]\right\}{\sigma }_{ui}^2, \\ U\left({\phi }_w\right)=\sum _{i=1}^n{x}_{3i}^{\top }\left\{\frac{1}{2{\sigma }_i^2}\left[\frac{e_i^2}{{\sigma }_i^2}+\frac{e_i{\lambda }_i}{{\sigma }_i}{A}_i\left(2+{\lambda }_i^2\right)-1\right]\right\}{\sigma }_{wi}^2+ \\ +\sum _{i=1}^n{\tilde{x}}_{3i}^{\top }\left\{\frac{1}{2}{\eta }^{\top }{\epsilon }_i\frac{{\sigma }_{wi}}{{\sigma }_{cw}}\left[\frac{e_i}{{\sigma }_i^2}+\frac{{\lambda }_i}{{\sigma }_i}{A}_i\right]\right\}, \end{gathered}$$ (9)

where $A_i=\frac{\varphi \left(a_i\right)}{\Phi \left(a_i\right)},a_i=-\frac{e_i{\lambda }_i}{{\sigma }_i}$, and ${\tilde{x}}_{3i}$ is x _3i , except for the null-intercept component.

From (8), when assuming an exponential distribution for u _i , the gradient is

$$\begin{gathered} U\left(\delta \right)=\sum _{i=1}^n{Z}_i^{\top }{\epsilon }_i{\Omega }^{-1}-\sum _{i=1}^n{Z}_i^{\top }\left\{\frac{{\sigma }_{wi}}{{\sigma }_{cw}}\left[\frac{1}{{\sigma }_{wi}}{B}_i-\frac{1}{{\sigma }_{ui}}\right]\right\}{\eta }^{\top }, \\ U\left(\eta \right)=\sum _{i=1}^n{\epsilon }_i^{\top }\left\{\frac{{\sigma }_{wi}}{{\sigma }_{cw}}\left[\frac{1}{{\sigma }_{wi}}{B}_i-\frac{1}{{\sigma }_{ui}}\right]\right\}, \\ U\left(\beta \right)=\sum _{i=1}^n{x}_{1i}^{\top }\left\{\frac{1}{{\sigma }_{wi}}{B}_i-\frac{1}{{\sigma }_{ui}}\right\}, \\ U\left({\phi }_u\right)=\sum _{i=1}^n{x}_{2i}^{\top }\left\{\frac{1}{2{\sigma }_{ui}^2}\left[\frac{{\sigma }_{wi}}{{\sigma }_{ui}}{B}_i-\frac{{\sigma }_{wi}^2}{{\sigma }_{ui}^2}-\frac{e_i}{{\sigma }_{ui}}-1\right]\right\}{\sigma }_{ui}^2, \\ U\left({\phi }_w\right)=\sum _{i=1}^n{x}_{3i}^{\top }\left\{\frac{1}{2{\sigma }_{ui}^2}+\frac{1}{2{\sigma }_{wi}}{B}_i\left[\frac{e_i}{{\sigma }_{wi}^2}-\frac{1}{{\sigma }_{ui}}\right]\right\}{\sigma }_{wi}^2+ \\ +\sum _{i=1}^n{\tilde{x}}_{3i}^{\top }\left\{\frac{1}{2}{\eta }^{\top }{\epsilon }_i\frac{{\sigma }_{wi}}{{\sigma }_{cw}}\left[\frac{1}{{\sigma }_{wi}}{B}_i-\frac{1}{{\sigma }_{ui}}\right]\right\}, \end{gathered}$$ (10)

where $B_i=\frac{\varphi \left(b_i\right)}{\Phi \left(b_i\right)},b_i=\frac{-e_i-{\sigma }_{wi}^2/{\sigma }_{ui}}{{\sigma }_{wi}}$, and ${\tilde{x}}_{3i}$ is x _3i , except for the null-intercept component.

From (8), when assuming a truncated normal distribution for u _i , the gradient is

$$\begin{gathered} U\left(\delta \right)=\sum _{i=1}^n{Z}_i^{\top }{\epsilon }_i{\Omega }^{-1}-\sum _{i=1}^n{Z}_i^{\top }\left\{\frac{{\sigma }_w}{{\sigma }_{cw}}\left[\frac{e_i+{\mu }_i}{{\sigma }^2}+\frac{\lambda }{\sigma }{D}_i\right]\right\}{\eta }^{\top }, \\ U\left(\eta \right)=\sum _{i=1}^n{\epsilon }_i^{\top }\left\{\frac{{\sigma }_w}{{\sigma }_{cw}}\left[\frac{e_i+{\mu }_i}{{\sigma }^2}+\frac{\lambda }{\sigma }{D}_i\right]\right\}, \\ U\left(\beta \right)=\sum _{i=1}^n{x}_{1i}^{\top }\left\{\frac{e_i+{\mu }_i}{{\sigma }^2}+\frac{\lambda }{\sigma }{D}_i\right\}, \\ U\left(\tau \right)=\sum _{i=1}^n{x}_{2i}^{\top }\left\{\frac{1}{\lambda \sigma }{D}_i-\frac{\sqrt{{\lambda }^{-2}+1}}{\sigma }{C}_i-\frac{e_i+{\mu }_i}{{\sigma }^2}\right\}, \\ U\left({\phi }_u\right)=\sum _{i=1}^n\left\{\frac{1}{2{\sigma }_u^2}\frac{{\mu }_i}{{\sigma }_u}{C}_i+\frac{1}{2{\sigma }^2}\left[\frac{(e_i+{\mu }_i{)}^2}{{\sigma }^2}-\frac{1}{\lambda \sigma }{D}_i\left(2{\mu }_i+e_i+\frac{{\mu }_i}{{\lambda }^2}\right)-1\right]\right\}{\sigma }_u^2, \\ U\left({\phi }_w\right)=\sum _{i=1}^n\left\{\frac{1}{2{\sigma }^2}\left[\frac{(e_i+{\mu }_i{)}^2}{{\sigma }^2}+\frac{\lambda }{\sigma }{D}_i\left({\mu }_i+2e_i+e_i{\lambda }^2\right)-1\right]\right\}{\sigma }_w^2, \end{gathered}$$ (11)

where $C_i=\frac{\varphi \left(c_i\right)}{\Phi \left(c_i\right)},c_i=\frac{{\mu }_i}{{\sigma }_u},D_i=\frac{\varphi \left(d_i\right)}{\Phi \left(d_i\right)}$, and $d_i=\frac{{\mu }_i}{\sigma \lambda }-\frac{e_i\lambda }{\sigma }$.

Alternatively, for computationally difficult cases, as suggested by Kutlu (2010); Karakaplan & Kutlu (2015); Amsler et al. (2016), it is possible to use a two-step maximum likelihood estimation method as in Murphy & Topel (2002). In this methodology, the parameter vector θ is estimated based on a two-step maximum likelihood estimation method (separately), characterizing a limited information maximum likelihood (LIML) approach.

Besides being easier to implement, the two-step estimation process can be extended to accommodate linear or non-linear regression models by instrumental variables. Thus, in the first stage, $\ln L_x\left(\theta \right)=\ln L_1\left({\theta }_1\right)$ is maximized in relation to its relevant parameters. In the second stage, conditional on the parameters estimated in the first stage, ln $\ln L_{y|x}\left(\theta \right)=\ln L_2\left({\theta }_2\right|{\theta }_1)$ is maximized.

The model of the second stage is

$y_i=x_{1i}^{\top }\beta +{\eta }^{\top }{\hat{\epsilon }}_i+e_i,$ (12)

where $e_i=w_i-u_i$ and ${\hat{\epsilon }}_i$ are the estimates of the first stage residuals obtained from the equation ${\hat{\epsilon }}_i=x_i-Z_i\hat{\delta }$ using ordinary least squares. Moreover, we can test the coefficients of the terms ${\hat{\epsilon }}_i$ for the presence of endogeneity by testing the null hypothesis that η = 0. In this structure, a variable is endogenous if it is not independent of v _i .

Therefore, $e_i=y_i-x_{1i}^{\top }\beta -{\eta }^{\top }{\hat{\epsilon }}_i$ and the other components are expressed in (8). As in the model proposed by Karakaplan & Kutlu (2015), in this approach, x follows a multivariate normal distribution if the number of endogenous variables is greater than one, and univariate normal otherwise. Whereas y|x is specified as normal/half-normal, normal/exponential or normal/truncated normal distribution.

A disadvantage compared to the one-step procedure is that although the two-step estimation leads to consistent estimation of θ ₂, the variance-covariance matrix estimated for y|x needs adjust. Due to the variability in ${\hat{\theta }}_1$, since ${\hat{\theta }}_1$ is an estimate of θ ₁ rather than its actual value. However, this approach presents fewer convergence problems.

Consequently, the two-step estimator provides incorrect and inconsistent standard errors, and a correction of these errors is required. To this end, an analytical approach is possible, as proposed by Murphy & Topel (2002). If the standard regularity conditions are met for both functions, then the two-step maximum likelihood estimator of θ ₂ is consistent and asymptotically normally distributed with a variance-covariance matrix

${\mathbf{V}}_{\mathbf{2}}^{\mathbf{*}}={\mathbf{V}}_{\mathbf{2}}\mathbf{+}{\mathbf{V}}_{\mathbf{2}}\mathbf{(}{\mathbf{CV}}_{\mathbf{1}}{\mathbf{C}}^{\mathbf{\top }}\mathbf{-}{\mathbf{RV}}_{\mathbf{1}}{\mathbf{C}}^{\mathbf{\top }}\mathbf{-}{\mathbf{CV}}_{\mathbf{1}}{\mathbf{R}}^{\mathbf{\top }}\mathbf{)}{\mathbf{V}}_{\mathbf{2}},$ (13)

where

$$\begin{gathered} {\mathbf{V}}_{\mathbf{1}}=(q\times q)\text{Asymptotic variance matrix of}{\hat{\theta }}_1\text{based on}\ln L_1\left({\theta }_1\right), \\ {\mathbf{V}}_{\mathbf{2}}=(p\times p)\text{Asymptotic variance matrix of}{\hat{\theta }}_2\text{based on}\ln L_2\left({\theta }_2\right|{\theta }_1), \\ \mathbf{C}=(p\times q)\text{matrix given by}E\left[\left(\frac{\partial \ln L_2}{\partial {\theta }_2}\right)\left(\frac{\partial \ln L_2}{\partial {\theta }_1^{\top }}\right)\right], \\ \mathbf{R}=(p\times q)\text{matrix given by}E\left[\left(\frac{\partial \ln L_2}{\partial {\theta }_2}\right)\left(\frac{\partial \ln L_1}{\partial {\theta }_1^{\top }}\right)\right]. \end{gathered}$$ (14)

The matrices V ₁ and V ₂ are estimated by the respective uncorrected variance-covariance matrices, typically by the inverse matrices of negative second derivatives. At the same time, the matrices C and R are estimated by summing the individual observations on the cross products of the derivatives.

A log-likelihood is assumed to exist for the first model, ln L ₁(θ1), as well as a conditional loglikelihood for the second (primary) model of interest, namely ln L ₂(θ ₂ |θ ₁). The component matrices of the Murphy-Topel estimator are estimated by the evaluation of the formulas in the maximum likelihood estimates for ${\hat{\theta }}_1$ and ${\hat{\theta }}_2$. As such,

${\hat{\mathbf{V}}}_{\mathbf{2}}^{\mathbf{*}}=\frac{1}{n}\left[{\hat{\mathbf{V}}}_{\mathbf{2}}\mathbf{+}{\hat{\mathbf{V}}}_{\mathbf{2}}\mathbf{(}\hat{\mathbf{C}}{\hat{\mathbf{V}}}_{\mathbf{1}}{\hat{\mathbf{C}}}^{\mathbf{\top }}\mathbf{-}\hat{\mathbf{R}}{\hat{\mathbf{V}}}_{\mathbf{1}}{\hat{\mathbf{C}}}^{\mathbf{\top }}\mathbf{-}\hat{\mathbf{C}}{\hat{\mathbf{V}}}_{\mathbf{1}}{\hat{\mathbf{R}}}^{\mathbf{\top }}\mathbf{)}{\hat{\mathbf{V}}}_{\mathbf{2}}\right],$ (15)

wherein

$$\begin{gathered} {\hat{\mathbf{V}}}_{\mathbf{1}}={\left[-\frac{1}{n}\sum _{i=1}^n\left(\frac{{\partial }^2\ln f_{1i}}{\partial {\hat{\theta }}_1\partial {\hat{\theta }}_1^{\top }}\right)\right]}^{-1},\hat{\mathbf{C}}=\frac{1}{n}\sum _{i=1}^n\left(\frac{\partial \ln f_{2i}}{\partial {\hat{\theta }}_2}\right)\left(\frac{\partial \ln f_{2i}}{\partial {\hat{\theta }}_1^{\top }}\right), \\ {\hat{\mathbf{V}}}_{\mathbf{2}}={\left[-\frac{1}{n}\sum _{i=1}^n\left(\frac{{\partial }^2\ln f_{2i}}{\partial {\hat{\theta }}_2\partial {\hat{\theta }}_2^{\top }}\right)\right]}^{-1},\hat{\mathbf{R}}=\frac{1}{n}\sum _{i=1}^n\left(\frac{\partial \ln f_{2i}}{\partial {\hat{\theta }}_2}\right)\left(\frac{\partial \ln f_{1i}}{\partial {\hat{\theta }}_1^{\top }}\right). \end{gathered}$$ (16)

By assuming exponential link function for ${\sigma }_{ui}^2$ and ${\sigma }_{wi}^2$, and ${\mu }_i\in \mathrm{ℝ}$, we obtain the following gradients:

The gradient of the two-step model, when assuming a half-normal distribution for u _i , is

$$\begin{gathered} U_1\left(\delta \right)=-2\sum _{i=1}^n{Z}_i^{\top }{\epsilon }_i, \\ {U}_2\left(\delta \right)=-\sum _{i=1}^n{Z}_i^{\top }\left\{\frac{e_i}{{\sigma }_i^2}+\frac{{\lambda }_i}{{\sigma }_i}{A}_i\right\}{\eta }^{\top }, \\ {U}_2\left(\eta \right)=\sum _{i=1}^n{\hat{\epsilon }}_i^{\top }\left\{\frac{e_i}{{\sigma }_i^2}+\frac{{\lambda }_i}{{\sigma }_i}{A}_i\right\}, \\ {U}_2\left(\beta \right)=\sum _{i=1}^n{x}_{1i}^{\top }\left\{\frac{e_i}{{\sigma }_i^2}+\frac{{\lambda }_i}{{\sigma }_i}{A}_i\right\}, \\ {U}_2\left({\phi }_u\right)=\sum _{i=1}^n{x}_{2i}^{\top }\left\{\frac{1}{2{\sigma }_i^2}\left[\frac{e_i^2}{{\sigma }_i^2}-\frac{e_i}{{\lambda }_i{\sigma }_i}{A}_i-1\right]\right\}{\sigma }_{ui}^2, \\ {U}_2\left({\phi }_w\right)=\sum _{i=1}^n{x}_{3i}^{\top }\left\{\frac{1}{2{\sigma }_i^2}\left[\frac{e_i^2}{{\sigma }_i^2}+\frac{e_i{\lambda }_i}{{\sigma }_i}(2+{\lambda }_i^2)A_i-1\right]\right\}{\sigma }_{wi}^2, \end{gathered}$$ (17)

where $A_i=\frac{\varphi \left(a_i\right)}{\Phi \left(a_i\right)}$ and $a_i=-\frac{e_i{\lambda }_i}{{\sigma }_i}$.

The gradient of the two-step model assuming an exponential distribution for u _i is

$$\begin{gathered} U_1\left(\delta \right)=-2\sum _{i=1}^n{Z}_i^{\top }{\epsilon }_i, \\ {U}_2\left(\delta \right)=-\sum _{i=1}^n{Z}_i^{\top }\left\{\frac{1}{{\sigma }_{wi}}{B}_i-\frac{1}{{\sigma }_{ui}}\right\}{\eta }^{\top }, \\ {U}_2\left(\eta \right)=\sum _{i=1}^n{\hat{\epsilon }}_i^{\top }\left\{\frac{1}{{\sigma }_{wi}}{B}_i-\frac{1}{{\sigma }_{ui}}\right\}, \\ {U}_2\left(\beta \right)=\sum _{i=1}^n{x}_{1i}^{\top }\left\{\frac{1}{{\sigma }_{wi}}{B}_i-\frac{1}{{\sigma }_{ui}}\right\}, \\ {U}_2\left({\phi }_u\right)=\sum _{i=1}^n{x}_{2i}^{\top }\left\{\frac{1}{2{\sigma }_{ui}^2}\left[\frac{{\sigma }_{wi}}{{\sigma }_{ui}}{B}_i-\frac{{\sigma }_{wi}^2}{{\sigma }_{ui}^2}-\frac{e_i}{{\sigma }_{ui}}-1\right]\right\}{\sigma }_{ui}^2, \\ {U}_2\left({\phi }_w\right)=\sum _{i=1}^n{x}_{3i}^{\top }\left\{\frac{1}{2{\sigma }_{ui}^2}+\frac{1}{2{\sigma }_{wi}}\left[\frac{e_i}{{\sigma }_{wi}^2}-\frac{1}{{\sigma }_{ui}}\right]B_i\right\}{\sigma }_{wi}^2, \end{gathered}$$ (18)

where $B_i=\frac{\varphi \left(b_i\right)}{\Phi \left(b_i\right)}$ and $b_i=\frac{-e_i-{\sigma }_{wi}^2/{\sigma }_{ui}}{{\sigma }_{wi}}$.

The gradient of the two-step model assuming a truncated normal distribution for u _i is

$$\begin{gathered} U_1\left(\delta \right)=-2\sum _{i=1}^n{Z}_i^{\top }{\epsilon }_i, \\ {U}_2\left(\delta \right)=-\sum _{i=1}^n{Z}_i^{\top }\left\{\frac{e_i+{\mu }_i}{{\sigma }^2}+\frac{\lambda }{\sigma }{D}_i\right\}{\eta }^{\top }, \\ {U}_2\left(\eta \right)=\sum _{i=1}^n{\hat{\epsilon }}_i^{\top }\left\{\frac{e_i+{\mu }_i}{{\sigma }^2}+\frac{\lambda }{\sigma }{D}_i\right\}, \\ {U}_2\left(\beta \right)=\sum _{i=1}^n{x}_{1i}^{\top }\left\{\frac{e_i+{\mu }_i}{{\sigma }^2}+\frac{\lambda }{\sigma }{D}_i\right\}, \\ {U}_2\left(\tau \right)=\sum _{i=1}^n{x}_{2i}^{\top }\left\{\frac{1}{\lambda \sigma }{D}_i-\frac{\sqrt{{\lambda }^{-2}+1}}{\sigma }{C}_i-\frac{e_i+{\mu }_i}{{\sigma }^2}\right\}, \\ {U}_2\left({\phi }_u\right)=\sum _{i=1}^n\left\{\frac{1}{2{\sigma }_u^2}\frac{{\mu }_i}{{\sigma }_u}{C}_i+\frac{1}{2{\sigma }^2}\left[\frac{(e_i+{\mu }_i{)}^2}{{\sigma }^2}-\frac{1}{\lambda \sigma }{D}_i\left(2{\mu }_i+e_i+\frac{{\mu }_i}{{\lambda }^2}\right)-1\right]\right\}{\sigma }_u^2, \\ {U}_2\left({\phi }_w\right)=\sum _{i=1}^n\left\{\frac{1}{2{\sigma }^2}\left[\frac{(e_i+{\mu }_i{)}^2}{{\sigma }^2}+\frac{\lambda }{\sigma }{D}_i\left({\mu }_i+2e_i+e_i{\lambda }^2\right)-1\right]\right\}{\sigma }_w^2, \end{gathered}$$ (19)

were $C_i=\frac{\varphi \left(c_i\right)}{\Phi \left(c_i\right)},c_i=\frac{{\mu }_i}{{\sigma }_u},D_i=\frac{\varphi \left(d_i\right)}{\Phi \left(d_i\right)}$ and $d_i=\frac{{\mu }_i}{\sigma \lambda }-\frac{e_i\lambda }{\sigma }$.

After obtaining estimates of the model parameters by maximum likelihood estimation, the next step is to predict the technical efficiency of each producer, ℰ_i = exp(−u _i ). A natural predictor for that amount is ${\hat{\mathcal{E}}}_i=\exp (-{\hat{u}}_i)$. However, Battese & Coelli (1988) used f(u _i |e _i ) to derive an alternative predictor, which was modified by Battese & Coelli (1995) to take into account the heteroscedasticity that may exist regarding the error components. This alternative predictor is

${\hat{\mathcal{E}}}_i=E\left\{\exp \right(-u_i\left)\right|e_i\}=\left\{\frac{\Phi ({\mu }_i^{*}/{\sigma }_i^{*}-{\sigma }_i^{*})}{\Phi ({\mu }_i^{*}/{\sigma }_i^{*})}\exp \left(\frac{1}{2}{\sigma }_i^{*2}-{\mu }_i^{*}\right)\right\},$ (20)

where ${\mu }_i^{*}$ and ${\sigma }_i^{*}$ vary according to the specification of u _i .

For the normal/half-normal model, ${\mu }_i^{*}$ and ${\sigma }_i^{*}$ are

$$\begin{gathered} {\mu }_i^{*}=\frac{-e_i{\sigma }_{ui}^2}{{\sigma }_i^2}, \\ {\sigma }_i^{*}=\frac{{\sigma }_{wi}{\sigma }_{ui}}{{\sigma }_i}. \end{gathered}$$ (21)

For the normal/exponential model, ${\mu }_i^{*}$ and ${\sigma }_i^{*}$ are

$$\begin{gathered} {\mu }_i^{*}=-e_i-\frac{{\sigma }_{wi}^2}{{\sigma }_{ui}}, \\ {\sigma }_i^{*}={\sigma }_{wi}. \end{gathered}$$ (22)

For the normal/truncated normal model, ${\mu }_i^{*}$ and ${\sigma }_i^{*}$ are

$$\begin{gathered} {\mu }_i^{*}=\frac{-e_i{\sigma }_u^2+{\mu }_i{\sigma }_w^2}{{\sigma }^2}, \\ {\sigma }^{*}=\frac{{\sigma }_w{\sigma }_u}{\sigma }. \end{gathered}$$ (23)

Battese & Coelli (1988) argue that, since the production function is usually defined by the logarithm of the production (ln y_i ), the technical efficiency for the i-th producer should be defined as $E\left\{\exp \right(-u_i\left)\right|e_i\}$. This predictor is optimal in terms of minimizing the mean squared prediction error.

As the data are in log terms, ${\hat{\mathcal{E}}}_i=E\left\{\exp \right(-u_i\left)\right|e_i\}$ is a measure of the percentage by which a unit fails to reach the frontier - the ideal production rate. Thus, the closer to one ${\hat{\mathcal{E}}}_i$ is, the closer the producer is to achieving optimal production, with the technology incorporated into the production function.

5 RESULTS

Following the approaches present in Section 4, we fitted six models to the data described in Section 3. Thus, we estimated the model parameters in one or two-step - full information maximum likelihood (FIML) or limited information maximum likelihood (LIML) - assuming a half-normal, exponential, or truncated normal distribution for the inefficiency term, u _i . For that, we postulated a Cobb-Douglas representation in a typical stochastic frontier approach under the endogeneity assumption of technical assistance and credit access variables.

Table 1 shows some goodness of fit measures of the one and two-step models. More specifically, it presents the log-likelihood values (ln L), the Akaike information criterion (AIC), and the Bayesian information criterion (BIC) of the models, as well as the Pearson correlation, the bias, and the root mean square error (RMSE) values between the observed and estimated values of the response variable, ln(income).

Table 1
- Goodness of fit measures.

Approach	Distribution of u	ln L	AIC	BIC	Cor(y, $\hat{\mathbf{y}}$ )	Bias	RMSE
FIML	Half-normal	-874.8	1827.5	2081.4	0.8792	0.19	1.135
	Exponential	-1330.6	2739.3	2993.2	0.8769	0.32	1.138
	Truncated normal	-432.2	944.5	1204.9	0.7315	2.87	2.895
LIML	Half-normal	-892.1	1814.3	1911.9	0.8792	0.20	1.126
	Exponential	-1284.7	2599.3	2697.0	0.8772	0.26	1.199
	Truncated normal	-435.7	903.5	1007.6	0.7636	2.88	2.901

When considering the exponential distribution for u _i , we did not obtain convergence with the BFGS optimization method, which uses analytic gradients. Therefore, we used the Nelder-Mead method, which does not use analytical gradients and is therefore slower. Until convergence, the normal/exponential model estimates using the FIML approach required 53280 iterations, whereas the LIML approach required only 4514. In both cases, compared to the other two specifications for u _i , the convergence rate when assuming an exponential distribution is considerably lower.

As for the goodness of fit measures, the normal/truncated normal models presented the lowest AIC and BIC and highest log-likelihood values, thus being the models better fitted by these criteria. However, the normal/half-normal models fitted better when considering the criteria of greater Pearson correlation and lower bias and REQM values. Additionally, all six models presented Pearson correlations between the observed and estimated values of the assumed endogenous variables (techassist and financing) greater than 0.8. Hence, indicating suitable fit of this variables from the linear regressions by instrumental variables. Based on these results and on the principle of parsimony, we selected the normal/half-normal models with their respective parameters estimated at one or two-step to model the total rural gross income of the Brazilian municipalities.

Therefore, the normal/half-normal stochastic production frontier model fitted is

$$\begin{gathered} \ln \left(incom{e}_i\right)={\beta }_0+{\beta }_1\ln \left(lan{d}_i\right)+{\beta }_2\ln \left(labo{r}_i\right)+{\beta }_3\ln \left(techinput{s}_i\right)+{\beta }_{4}regio{n}_{nort{h}_i}+ \\ +{\beta }_{5}regio{n}_{northeas{t}_i}+{\beta }_{6}regio{n}_{southeas{t}_i}+{\beta }_{7}regio{n}_{sout{h}_i}+v_i-u_i, \\ \ln \left({\sigma }_{ui}^2\right)={\phi }_{u0}+{\phi }_{u1}techassis{t}_i+{\phi }_{u2}financin{g}_i+{\phi }_{u3}gin{i}_i, \\ \ln \left({\sigma }_{wi}^2\right)={\phi }_{w0}, \\ techassis{t}_i={\delta }_0+{\delta }_1\ln \left(lan{d}_i\right)+{\delta }_2\ln \left(labo{r}_i\right)+{\delta }_3\ln \left(techinput{s}_i\right)+{\delta }_{4}regio{n}_{nort{h}_i}+ \\ +{\delta }_{5}regio{n}_{northeas{t}_i}+{\delta }_{6}regio{n}_{southeas{t}_i}+{\delta }_{7}regio{n}_{sout{h}_i}+{\delta }_{8}socia{l}_i+ \\ +{\delta }_{9}demographi{c}_i+{\delta }_{10}environmen{t}_i+{\delta }_{11}gin{i}_i+{\epsilon }_{1i}, \\ financin{g}_i={\gamma }_0+{\gamma }_1\ln \left(lan{d}_i\right)+{\gamma }_2\ln \left(labo{r}_i\right)+{\gamma }_3\ln \left(techinput{s}_i\right)+{\gamma }_{4}regio{n}_{nort{h}_i}+ \\ +{\gamma }_{5}regio{n}_{northeas{t}_i}+{\gamma }_{6}regio{n}_{southeas{t}_i}+{\gamma }_{7}regio{n}_{sout{h}_i}+{\gamma }_{8}socia{l}_i+ \\ +{\gamma }_{9}demographi{c}_i+{\gamma }_{10}environmen{t}_i+{\gamma }_{11}gin{i}_i+{\epsilon }_{2i}, \end{gathered}$$

where the Center-West region is the base level.

Tables 2 and 3 provide estimates of the normal/half-normal full information and limited information models, respectively. In the ^Appendix (Tables 6 and 7) are the results from the linear regressions by instrumental variables of these models.

Table 2
- Full information maximum likelihood estimation of the normal/half-normal model.

Variable	Coefficient	Standard error	z	P-value	Lower limit	Upper limit
Frontier
constant	2.477	0.125	19.852	0.000	2.233	2.722
ln(land)	0.110	0.014	7.639	0.000	0.082	0.139
ln(labor)	0.206	0.011	17.945	0.000	0.184	0.229
ln(techinputs)	0.509	0.018	28.977	0.000	0.475	0.544
regionnorth	0.066	0.055	1.201	0.230	-0.042	0.174
regionnortheast	0.129	0.053	2.425	0.015	0.025	0.233
regionsoutheast	0.242	0.046	5.211	0.000	0.151	0.333
regionsouth	0.335	0.048	7.003	0.000	0.241	0.429
ln( ${\mathit{\sigma }}_{\mathbf{u}}^{\mathbf{2}}$ )
constant	5.637	0.201	28.040	0.000	5.243	6.031
techassist	0.931	0.525	1.772	0.076	-0.099	1.961
financing	-2.931	0.579	-5.061	0.000	-4.067	-1.796
gini	-9.938	0.299	-33.267	0.000	-10.523	-9.352
ln( ${\mathit{\sigma }}_{\mathbf{w}}^{\mathbf{2}}$ )
constant	-1.078	0.023	-47.806	0.000	-1.122	-1.034
η techassist
constant	0.890	0.079	11.194	0.000	0.734	1.045
η financing
constant	0.200	0.081	2.484	0.013	0.042	0.358
${\mathit{\sigma }}_{\mathbf{w}}^{\mathbf{2}}$	0.340	0.008	44.359	0.000	0.325	0.355

Table 3
- Limited information maximum likelihood estimation of the normal/half-normal model.

Variable	Coefficient	Standard error	z	P-value	Lower limit	Upper limit
Frontier
constant	2.552	0.126	20.306	0.000	2.306	2.799
ln(land)	0.111	0.015	7.521	0.000	0.082	0.140
ln(labor)	0.202	0.012	16.813	0.000	0.178	0.226
ln(techinputs)	0.507	0.018	28.133	0.000	0.472	0.543
regionnorth	0.062	0.055	1.123	0.261	-0.046	0.170
regionnortheast	0.131	0.053	2.455	0.014	0.026	0.235
regionsoutheast	0.231	0.046	4.995	0.000	0.140	0.321
regionsouth	0.319	0.047	6.736	0.000	0.226	0.412
ln( ${\mathit{\sigma }}_{\mathbf{u}}^{\mathbf{2}}$ )
constant	6.337	0.594	10.662	0.000	5.172	7.502
techassist	-0.182	0.578	-0.315	0.753	-1.314	0.950
financing	-2.215	0.617	-3.592	0.000	-3.423	-1.006
gini	-10.441	0.841	-12.420	0.000	-12.089	-8.794
ln( ${\mathit{\sigma }}_{\mathbf{w}}^{\mathbf{2}}$ )
constant	-1.065	0.023	-45.429	0.000	-1.111	-1.019
η techassist
constant	0.664	0.080	8.335	0.000	0.508	0.820
η financing
constant	0.238	0.080	2.986	0.003	0.082	0.395
${\mathit{\sigma }}_{\mathbf{w}}^{\mathbf{2}}$	0.345	0.008	43.938	0.000	0.329	0.360

Note that, as expected for this type of model, the components of expenditure on land, labor, and capital have significant positive effects on income (Tables 2 and 3). In addition, credit access (financing) and income concentration (gini) have significant negative effects on the ${\sigma }_{ui}^2$ function for the technical inefficiency level of agricultural properties. Indicating that greater access to rural credit and income concentration reduces the inefficiency of agricultural properties. In contrast, technical assistance (techassist) does not have a significant effect at 5%. This result is due to market imperfections, here represented by income concentration, which prevents technical assistance from being significant. Note the evidence of endogeneity in both cases (${\hat{\eta }}_{financing}$ and ${\hat{\eta }}_{techassist}$ with p-value < 0.05).

It is noteworthy that Souza & Gomes (2018) achieved significance for technical assistance by removing income concentration from the analysis. They concluded that the social indicator is the key variable to reduce inefficiency and reported technical assistance as a significant part of rural extension positively affecting income. In addition, improving the social indicator will facilitate access to technical assistance, thus creating a positive synergistic effect on income, reducing income concentration.

In the present work, we found that the technical assistance indicator is relatively too low for the Northern and Northeastern regions - in general, the values are less than half of the corresponding values for the other regions. Thus, public policies should be oriented to improve this indicator, especially in these regions.

Table 4 shows the results of the Wald and likelihood ratio tests on the presence of endogeneity. Note that we rejected the null hypothesis of exogeneity in both approaches - evidence of endogeneity. Therefore, to obtain consistent parameter estimates, we need fit models that take endogeneity into account.

Table 4
- Wald and likelihood ratio tests for endogeneity.

Approach	Wald	Likelihood ratio	P-value
FIML	201	342	0.000
LIML	149	208	0.000

Table 5 summarizes the relative importance of the production factors, including returns to scale for the one and two-step models. In both cases, we get decreasing returns to scale. Furthermore, capital (technological inputs) dominates the production function, followed by labor and land, showing that capital as input has a greater influence on production, corroborating the literature.

Table 5
- Relative elasticities and returns to scale.

Production factor	Coefficient	Standard error	t	P-value	Lower limit	Upper limit
FIML
Labor	0.25	0.01	18.30	0.00	0.22	0.28
Land	0.13	0.02	7.69	0.00	0.10	0.17
Capital	0.62	0.02	31.01	0.00	0.58	0.66
Returns to scale	0.83	0.01
LIML
Labor	0.25	0.01	18.13	0.00	0.22	0.27
Land	0.14	0.02	7.91	0.00	0.10	0.17
Capital	0.62	0.02	31.47	0.00	0.58	0.66
Returns to scale	0.82	0.01

Figure 1 illustrates the box plots for the normalized classifications of the technical efficiency measurements by region (${\hat{\mathcal{E}}}_i$) predicted by the normal/half-normal models using FIML and LIML approaches. We have that efficiency differs significantly by region. Note the predominance of the Center-West region over the others, followed by the Southeast and South, and that the North and Northeast have the lowest efficiency levels.

Figure 1
- Box plots for the technical efficiency measurements by region predicted by the normal/half-normal models using FIML and LIML approaches.

Estimates of the one and two-step procedures (FIML or LIML approaches) and predicted technical efficiencies are similar when using a normal/half-normal specification. Consequently, under standard regularity conditions, we recommend using the one-step procedure instead of the twostep one. In such conditions, the FIML estimator is more efficient than the LIML estimator and generally produces the lowest standard deviations. In contrast, we recommend using the twostep procedure for computationally intensive cases or when the one-step procedure reaches no convergence.

6 CONCLUSIONS

This paper implements the prediction of the producers’ technical efficiency level from stochastic production frontier models with endogenous and exogenous variables and heteroscedastic error terms through one and two-step maximum likelihood estimation, based on Karakaplan & Kutlu (2015). We consider three main specifications for the inefficiency term (half-normal, exponential, and truncated normal distributions). We also derived the analytic gradients of these models, which are not yet available in the literature and can provide better performance at a reasonable computational time cost. Moreover, we implemented functions in the R language to the methodology presented here.

Additionally, we apply the models to municipal data from the Brazilian agricultural census. The application favored the use of the proposed regression models. The results from the normal/halfnormal stochastic production frontier models under endogeneity are remarkably similar. Therefore, if there is convergence in the one and two-step models, then the one-step maximum likelihood estimator is recommended due to its better efficiency - smaller standard errors.

It is important to note that the correction of the variance-covariance matrix made in the two-step estimation method can change the significance of important variables compared to those in the one-step. Thus, it can change the expected technical efficiencies, especially when applying a normal/exponential model, which usually presents convergence problems.

The production function estimation is dominated by capital (technological inputs), followed by labor and land. Production shows decreasing returns to scale. Credit access and technical assistance are endogenous, and income concentration seems to impede productive inclusion through the more intensive use of technology.

Further studies are required to implement functions in R allowing: alternative parameterizations for the inefficiency term, such as the gamma distribution; nonlinear regressions by instrumental variables for the assumed endogenous variables; and different diagnostic analyses, such as residuals analysis. However, at the moment, a routine in the R language is available for the approaches described in this paper.

References

• AIGNER D, LOVELL CK & SCHMIDT P. 1977. Formulation and estimation of stochastic frontier production function models. Journal of econometrics, 6(1): 21-37.
• AMSLER C, PROKHOROV A & SCHMIDT P. 2016. Endogeneity in stochastic frontier models. Journal of Econometrics, 190(2): 280-288.
• ANDRADE BB & SOUZA GS. 2017. Likelihood computation in the normal-gamma stochastic frontier model. Computational Statistics, pp. 1-16.
• ANDRADE BB & SOUZA GS. 2019. The EM algorithm for standard stochastic frontier models. Pesquisa Operacional, 39: 361-378.
• BATTESE GE & COELLI TJ. 1988. Prediction of firm-level technical efficiencies with a generalized frontier production function and panel data. Journal of econometrics , 38(3): 387-399.
• BATTESE GE & COELLI TJ. 1992. Frontier production functions, technical efficiency and panel data: with application to paddy farmers in India. Journal of productivity analysis, 3(1-2): 153- 169.
• BATTESE GE & COELLI TJ. 1995. A model for technical inefficiency effects in a stochastic frontier production function for panel data. Empirical economics, 20(2): 325-332.
• CAZALS C, FE` VE F, FLORENS JP & SIMAR L. 2016. Nonparametric instrumental variables estimation for efficiency frontier. Journal of econometrics , 190(2): 349-359.
• GREENE WH. 1990. A gamma-distributed stochastic frontier model. Journal of econometrics , 46(1-2): 141-163.
• GRIFFITHS WE & HAJARGASHT G. 2016. Some models for stochastic frontiers with endogeneity. Journal of Econometrics , 190(2): 341-348.
• GUAN Z, KUMBHAKAR SC, MYERS RJ & LANSINK AO. 2009. Measuring excess capital capacity in agricultural production. American Journal of Agricultural Economics, 91(3): 765- 776.
• KARAKAPLAN MU. 2017. Fitting endogenous stochastic frontier models in Stata. The Stata Journal, 17(1): 39-55.
• KARAKAPLAN MU & KUTLU L. 2015. Handling endogeneity in stochastic frontier analysis. Available at SSRN 2607276, .
• KUMBHAKAR SC & LOVELL CK. 2003. Stochastic frontier analysis. Cambridge university press.
• KUMBHAKAR SC, PARMETER CF & ZELENYUK V. 2020. Stochastic frontier analysis : Foundations and advances I. Handbook of production economics, pp. 1-40.
• KUTLU L. 2010. Battese-Coelli estimator with endogenous regressors. Economics Letters, 109(2): 79-81.
• MEEUSEN W & VAN DEN BROECK J. 1977. Efficiency estimation from Cobb-Douglas production functions with composed error. International economic review, pp. 435-444.
• MURPHY KM & TOPEL RH. 2002. Estimation and inference in two-step econometric models. Journal of Business & Economic Statistics, 20(1): 88-97.
• MUTTER RL, GREENE WH, SPECTOR W, ROSKO MD & MUKAMEL DB. 2013. Investigating the impact of endogeneity on inefficiency estimates in the application of stochastic frontier analysis to nursing homes. Journal of productivity analysis , 39(2): 101-110.
• PROKHOROV A, TRAN KC & TSIONAS MG. 2020. Estimation of semi- and nonparametric stochastic frontier models with endogenous regressors. Empirical Economics, 60(6): 3043-3068.
• SOUZA GS & GOMES EG. 2018. A stochastic production frontier analysis of the Brazilian agriculture in the presence of an endogenous covariate. In: International Conference on Operations Research and Enterprise Systems. pp. 3-14. Springer.
• SOUZA GS, GOMES EG & ALVES ERA. 2017. Conditional FDH efficiency to assess performance factors for Brazilian agriculture. Pesquisa Operacional , 37: 93-106.
• TRAN KC & TSIONAS EG. 2013. GMM estimation of stochastic frontier model with endogenous regressors. Economics Letters , 118(1): 233-236.
• TRAN KC & TSIONAS EG. 2015. Endogeneity in stochastic frontier models: Copula approach without external instruments. Economics Letters , 133: 85-88.
• TSIONAS M, IZZELDIN M, HENNINGSEN A & PARAVALOS E. 2021. Addressing endogeneity when estimating stochastic ray production frontiers: a Bayesian approach. Empirical Economics , pp. 1-19.

APPENDIX

Table 6
Instrumental variables regression of the one-step model.

Variable	Coefficient	Standard error	z	P-value	Lower limit	Upper limit
techassist
constant	-0.260	0.036	-7.277	0.000	-0.331	-0.190
ln(land)	0.005	0.004	1.247	0.212	-0.003	0.013
ln(labor)	0.010	0.003	3.071	0.002	0.004	0.017
ln(techinputs)	0.082	0.005	17.294	0.000	0.073	0.091
regionnorth	0.048	0.015	3.129	0.002	0.018	0.078
regionnortheast	0.048	0.015	3.201	0.001	0.019	0.078
regionsoutheast	0.071	0.013	5.292	0.000	0.045	0.097
regionsouth	0.165	0.014	11.542	0.000	0.137	0.194
social	0.406	0.023	17.754	0.000	0.361	0.451
demographic	-0.016	0.028	-0.560	0.575	-0.071	0.040
environment	0.040	0.034	1.180	0.238	-0.027	0.107
gini	-0.579	0.031	-18.884	0.000	-0.639	-0.519
financing
constant	-0.508	0.036	-13.939	0.000	-0.580	-0.437
ln(land)	0.027	0.004	6.831	0.000	0.019	0.035
ln(labor)	-0.005	0.003	-1.381	0.167	-0.011	0.002
ln(techinputs)	0.129	0.005	26.873	0.000	0.120	0.138
regionnorth	-0.076	0.015	-4.884	0.000	-0.106	-0.045
regionnortheast	-0.080	0.015	-5.241	0.000	-0.110	-0.050
regionsoutheast	-0.057	0.014	-4.223	0.000	-0.084	-0.031
regionsouth	0.104	0.015	7.168	0.000	0.076	0.133
social	0.156	0.024	6.590	0.000	0.109	0.202
demographic	-0.222	0.029	-7.576	0.000	-0.279	-0.165
environment	-0.398	0.036	-11.173	0.000	-0.467	-0.328
gini	-0.197	0.032	-6.087	0.000	-0.260	-0.133

Table 7
Instrumental variables regression of the two-step model.

Variable	Coefficient	Standard error	z	P-value	Lower limit	Upper limit
techassist
constant	-0.293	0.036	-8.148	0.000	-0.363	-0.222
ln(land)	0.003	0.004	0.864	0.388	-0.004	0.011
ln(labor)	0.004	0.003	1.167	0.243	-0.003	0.010
ln(techinputs)	0.079	0.005	16.716	0.000	0.070	0.088
regionnorth	0.055	0.015	3.598	0.000	0.025	0.085
regionnortheast	0.047	0.015	3.124	0.002	0.018	0.077
regionsoutheast	0.059	0.013	4.419	0.000	0.033	0.085
regionsouth	0.155	0.014	10.819	0.000	0.127	0.183
social	0.487	0.022	21.740	0.000	0.443	0.531
demographic	-0.003	0.029	-0.109	0.913	-0.060	0.054
environment	-0.018	0.035	-0.510	0.610	-0.086	0.051
gini	-0.426	0.029	-14.643	0.000	-0.483	-0.369
financing
constant	-0.521	0.036	-14.313	0.000	-0.593	-0.450
ln(land)	0.027	0.004	6.682	0.000	0.019	0.035
ln(labor)	-0.007	0.003	-2.186	0.029	-0.014	-0.001
ln(techinputs)	0.128	0.005	26.669	0.000	0.118	0.137
regionnorth	-0.073	0.015	-4.698	0.000	-0.103	-0.042
regionnortheast	-0.081	0.015	-5.269	0.000	-0.111	-0.051
regionsoutheast	-0.062	0.014	-4.588	0.000	-0.089	-0.036
regionsouth	0.100	0.015	6.879	0.000	0.072	0.129
social	0.189	0.023	8.320	0.000	0.145	0.234
demographic	-0.217	0.029	-7.359	0.000	-0.275	-0.159
environment	-0.421	0.035	-11.912	0.000	-0.491	-0.352
gini	-0.134	0.030	-4.527	0.000	-0.192	-0.076

Publication Dates

• Publication in this collection
  02 May 2022
• Date of issue
  2022

History

• Received
  11 Sept 2020
• Accepted
  18 Jan 2022