Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
79,429 characters · 8 sections · 65 citation commands
Maximum Likelihood Estimation of Stochastic Frontier Models with Endogeneity
\onehalfspacing
Endogeneity in the stochastic frontier framework has received increasing attention in recent work (kutlu2010,tran2013,tran2015,karakaplan2017,amsler2016,Lai2018, and kumbhakar2020a,kumbhakar2020b, for a review). Most contributions focus on the correlation between the regressors and the two-sided error component while ignoring the potential dependence between regressors and the stochastic inefficiency component. However, if producers have some information about their inefficiency level, they can use it to guide their choice of inputs and environmental variables (e.g., managerial characteristics). That is, there may be factors, observable to the firm but unobservable to the econometrician, which affect both the choice of regressors, and the level of inefficiency cazals2016.
In this paper, we consider a stochastic frontier model in which both the two-sided error and the stochastic inefficiency terms are allowed to be correlated with inputs and environmental variables. These endogenous variables are restricted to be continuous. The production frontier can be linear or nonlinear, and the inefficiency term satisfies the scaling property. That is, it can be decomposed into a stochastic efficiency term and a scaling function that depends on environmental variables alvarez2006. We achieve identification by allowing for a vector of control functions that fully captures the dependence between the composite error term and the endogenous variables.
To the best of our knowledge, models that explicitly allow for dependence between the stochastic inefficiency term, inputs, and environmental variables have only been studied by amsler2017.\footnote{karakaplan2017 do not directly consider the potential endogeneity of the inefficiency term. They instead model the potential dependence between the two-sided error term and the inefficiency term through observables, while we model such dependence through unobservables.} In their paper, the marginal distribution of the statistical noise is taken to be a normal distribution, and the marginal distribution of the stochastic inefficiency term to be a half-normal distribution. The two stochastic terms are potentially correlated. The dependence between observables and unobservables is modeled using copula functions, which are cleverly constructed from the marginal distributions of the unobservables. However, the likelihood function cannot be written in closed form, and the authors need to resort to simulations to obtain an estimator of the model's parameters. This approach prevents a clear analysis of identification, estimation, and inference. Moreover, simulated methods can be biased and have a higher variance in finite samples, especially when the number of simulations is not chosen appropriately with the sample size gourieroux1996. Finally, when both inputs and environmental variables are potentially correlated with the inefficiency term, they cannot obtain an estimator of technical efficiency.
Our framework seeks to avoid these potential pitfalls. In particular, we are able to obtain the maximum likelihood function in closed form. This allows us to further the analysis of identification of this model and propose a simple and computationally fast estimation of the model's parameters. We also offer a generalization of the battese1988 estimator of technical efficiency.
While our main statistical model is similar to the one in amsler2017, two fundamental assumptions deviate from their framework. First, we assume that the two-sided error term and the inefficiency term are independent conditional on a vector of control functions. This first assumption allows us to write the conditional density of the composite error as the convolution of the conditional densities of the statistical noise and the inefficiency term, respectively. Second, we assume that the conditional distribution of the baseline stochastic inefficiency term given the control functions is a folded normal distribution leone1961,sundberg1974. The latter assumption is convenient for two main reasons. On the one hand, it allows us to capture the dependence between inefficiency and endogenous variables through a vector of what we refer to as dependence parameters, $\rho_{U}$, which take values in the hypercube $[-1,1]$. These parameters measure the fraction of inefficiency observed by the producer but not by the econometrician, which may influence the choice of inputs and environmental factors and is confounded with the observed level of the regressors. On the other hand, the conditional normal-folded normal model provides a natural generalization of the normal half-normal model to the case when inputs and environmental variables are endogenous. This is because the folded normal pdf collapses to the density of a half-normal random variable when $\rho_{U} = \mathbf{0}$, where $\mathbf{0}$ is a vector of zeros. That is, when stochastic inefficiency is unobservable to producers and hence cannot influence their decision.
Our analysis of identification and estimation focuses on the dependence parameters $\rho_U$. Because of the properties of the folded normal distribution, only the magnitude of the components of $\rho_U$ is identified. However, their sign cannot be identified sundberg1974,schmidt1980.\footnote{The folded normal distribution can be thought of as a normal distribution “folded" at zero by taking the absolute value. Suppose we take a mean-zero normal random variable $\eta$, and then generate two standard normal random variables $U_1$ and $U_2$, which have correlation $-0.5$, and $0.5$ with $\eta$, respectively. When we “fold" both $U_1$ and $U_2$ by taking their absolute values, we have that $\vert U_1 \vert$ has the same conditional distribution of $\vert U_2 \vert$. Identification of the sign of $\rho_U$ is thus not feasible.} Hence, the likelihood function has two isolated maxima, which are symmetric about a local extremum at zero. We deal with this identification issue by imposing a sign normalization which amounts to restricting one of the components of $\rho_U$ to lay in the positive orthant. When $\rho_U = \mathbf{0}$, the likelihood function has a unique extremum. However, when $\rho_U = \mathbf{0}$ the score is identically equal to zero, and our model is not first-order identified. We are nonetheless able to show that our model is second-order identified.
Moreover, since one of the components of $\rho_U$ lays at the boundary of the parameters' space and the Hessian matrix is singular when $\rho_U = \mathbf{0}$, our estimator has a non-standard asymptotic distribution, and its rate of convergence is slower than $\sqrt{n}$, where $n$ is the sample size. We provide the asymptotic distribution of our estimator in all these cases following the framework of andrews1999 and rot2000. Finally, we briefly discuss potential ways to conduct inference on $\rho_U$.
Throughout the paper, we assume that parameters other than $\rho_U$ are first-order locally identified, and thus $\sqrt{n}$-estimable sargan1983. This assumption implies, in particular, that the variance of the inefficiency term is strictly positive. We defer to future research the study of this model when such an assumption fails lee1993.
The paper is structured as follows. In Section (ref), we discuss the statistical model and provide the main steps for the construction of the likelihood function. We further consider identification, estimation and inference. In Section (ref), we provide simulation evidence of the finite sample properties of our estimator. We show that our estimator performs better than the copula method of amsler2017 especially for estimating the variance of the stochastic inefficiency term. In Section (ref), we apply our methodology to the agricultural sector in Nepal. We show that accounting for endogeneity substantially changes the conclusions of the empirical analysis. In particular, our estimator detects considerable variation in the efficiency scores which is not found when regressors are taken to be exogenous.
We study a general version of the model usually considered in this literature. The logarithm of the output, $Y$, is determined by some known function, $m(\cdot,\cdot)$, which depends on a vector of $p \geq 1$ inputs, $X$, and parameters, $\beta$; and by a composite error term $\varepsilon = V - U$, where $V$ represents a stochastic component; and $U \geq 0$ is the so-called inefficiency term. We thus have
where $U$ captures the producer's shortfall from the production frontier.
Additionally, we fix $U = U_0g(Z,\delta)$, where $U_0\geq 0$ is a stochastic inefficiency component and $g(\cdot,\cdot)$ is a known strictly positive scaling function, which depends on some additional environmental variables $Z \in \mathbb{R}^k$, with $k \geq 0$, through parameters $\delta$ simar1994,alvarez2006. The scaling function further satisfies the normalization condition $g(0,\delta) = 1$. $X$ and $Z$ may have some common elements, but they must have at least one non-overlapping component.
Thus, we finally have
A maximum likelihood estimator of $(\beta,\delta)$ is based on the assumption that the composite error component $(V,U_0)$ is independent of $(X,Z)$, with $(U_0,V)$ mutually independent; $V$ following a normal distribution with a constant variance, and $U_0$ following a normal distribution truncated at $0$ aigner1977,schmidt1979,schmidt1980,horrace2005. While a consistent estimation of $(\beta,\delta)$ can also be obtained without these strong distributional assumptions simar1994,tran2013, these assumptions are necessary to learn something about the variance of the inefficiency term, $U_0$. We are often interested in estimating each producer's distance from the frontier battese1988. This can be easily done when the marginal distributions of $V$ and $U_0$ are taken to be known.
The literature has long recognized that inputs may be simultaneously chosen with the output, and thus potentially correlated with the composite error term mundlak1961,schmidt1984. Similarly, the producer may decide environmental variables depending on characteristics that are observable to her but not to the econometrician.
To deal with endogenous variables, we need a vector of instruments that are correlated with the endogenous components but independent of the composite error term amsler2016. To simplify our presentation, we take all variables in $(X,Z)$ to be endogenous. The extension to the case when we have some endogenous and some exogenous components can be handled similarly.
We consider the following auxiliary regression models
where $\eta = (\eta^\prime_X,\eta^\prime_Z)^\prime \in \mathbb{R}^{p+k}$ is a random vector of error components, and $W\in \mathbb{R}^q$ is a vector of instrumental variables, with $q \geq p + k$.
Our approach is based on a control function assumption. That is, we assume that all the dependence between $(X,Z)$ and $(V,U_0)$ is captured by $\eta$ newey1999,newey2009,wooldridge2009. Moreover, we assume that the instruments are strongly exogenous, that is, independent of the composite error term. Given a triplet of random variables $U_0$, $V$ and $\eta$, we use the notation $U_0 \upmodels V$ to indicate that $U_0$ is independent of $V$; and the notation $U_0 \upmodels V \vert \eta$ to indicate that $U_0$ is independent of $V$ conditional on $\eta$.
Our independence assumptions can be formally stated as follows:
Assumption (ref) implies strong exogeneity of the instruments; and that the control function, $\eta$, captures all the dependence between $(X,Z)$ and $(U_0,V)$. That is, $(X,Z) \upmodels (U_0,V) \vert \eta$.
Assumption (ref) implies that, if any dependence exists between $V$ and $U_0$, it has to happen through the vector $\eta$. This assumption reduces to the standard assumption of $U_0 \upmodels V$ when both $X$ and $Z$ are taken to be exogenous kumbhakar2003.
Assumptions (ref) and (ref) imply that \[ f_{V,U_0,\eta \vert W} (v,u,\eta \vert W) = f_{V,U_0,\eta} (v,u,\eta )= f_{V,\eta} (v ,\eta ) f_{U_0\vert \eta} (u \vert\eta ), \] where $f$ denotes a probability density function. To construct a maximum likelihood estimator (MLE), we let $\eta \sim N(0,\Sigma_\eta)$, where $\Sigma_\eta$ is a symmetric, positive definite covariance matrix. We also let $D_{\eta}$ be a $p+ k \times p + k$ diagonal matrix whose diagonal entries are the standard deviations of $\eta$. We can write that $\Sigma_\eta = D_{\eta}C_\eta D_{\eta}$, where $C_\eta$ is the symmetric, positive definite correlation matrix of the vector $\eta$. That is, a matrix with diagonal equal to $1$ and the other elements in the interval $(-1,1)$.
An additional requirement for the construction of a full information MLE is that \[
\sim N\left(
,
\right), \] where $\rho_{V}$ is a vector of correlation coefficients between $V$ and all components of $\eta$, and $\sigma^2_V$ is the variance of $V$ kutlu2010.
The main difficulty lies in the specification of the joint density of $(U_0,\eta)$ such that its marginal distributions are a half-normal and a joint normal, respectively, and the dependence between the two is captured by only one parameter. If one specifies a joint normal distribution for the random vector $(U^\ast_0,\eta)$ and then takes $U_0 = \vert U^\ast_0\vert$, the marginal distributions of $U_0$ and $\eta$ are the correct marginal distributions. This construction also creates dependence between $U_0$ and $\eta$. Here, we argue that the conditional distribution of $U_0$ given $\eta$ can be written in such a way that this dependence is captured by only one vector of parameters which, we refer to as dependence parameters, and we denote as $\rho_U$. Let \[
\sim N\left(
,
\right), \] where $\rho_{U}$ is a vector of correlations between $U^\ast_0$ and $\eta$, and $\sigma^2_U$ is the variance of $U^\ast_0$. The conditional density of $U^\ast_0$ given $\eta$ is \[ f_{U^\ast_0 \vert \eta } (u \vert\eta ) = \frac{1}{\sqrt{2\pi \sigma_U^2 (1 - \rho^\prime_{U}C^{-1}_\eta \rho_{U})}} \exp \left( -\frac{(u - \sigma_U \rho^\prime_{U} C^{-1}_\eta D^{-1}_\eta \eta)^2}{2 \sigma_U^2 (1 - \rho^\prime_{U}C^{-1}_\eta \rho_{U})}\right). \] Thus, we have that \[ P\left( U_0 \leq u \vert \eta \right) = P\left( U^\ast_0 \leq u \vert \eta \right) - P\left( U^\ast_0 \leq -u \vert \eta \right). \] Taking the derivative of the last equality with respect to $u$ on both sides, we obtain that the density of $U_0$ given $\eta$ is equal to \[ f_{U_0 \vert \eta }(u \vert \eta )= f_{U^\ast_0 \vert \eta } (u \vert\eta ) + f_{U^\ast_0 \vert \eta } (-u \vert\eta ). \] Therefore, the conditional density function of $U_0$ is
which is the pdf of a folded normal distribution leone1961. When we impose that $\rho_{U}$ is a vector of zeros, that is, when there is no dependence between the regressors and the inefficiency term, the conditional distribution in (ref) reduces to \[ f_{U_0}(u) = \frac{2}{\sqrt{2\pi\sigma^2_U}} \exp \left( -\frac{u^2}{2\sigma_U^2}\right), \] which is the density of a half-normal distribution. We also show in Appendix (ref) that the marginal density of $U_0$ obtained from this construction is a half-normal density.
Figure (ref) depicts the conditional folded normal pdf when $\eta$ is a bivariate random vector with unit variance and correlation coefficient equal to $0.5$, $\sigma^2_U = 2.752$, and $\rho_U = (0.5,0.5)^\prime$. For fixed parameters, the pdf is symmetric in $\eta$, in the sense that the shape of the density for $\eta =e$ is the same as for $\eta = -e$, for any real-valued vector $e$.
For a given $\eta$, this implies that the density is invariant to changes in sign of the vector of dependence parameters $\rho_U$. That is, the conditional density of $U_0$ generated under a certain dependence vector $\rho_U$ is equal to the conditional density of $U_0$ when the dependence vector is $-\rho_U$. This is a well-known equivalence property of the folded normal distribution sundberg1974.
The construction of the likelihood function is thus based on the following
Finally, because of Assumption (ref) and the strict positivity of the function $g(\cdot,\cdot)$, the conditional distribution of $U = U_0 g(Z,\delta)$ given $\eta$ can be written as
and it is therefore a simple scaled version of the distribution of $U_0$ given $\eta$, as in the exogenous case.
We follow the literature on stochastic frontier and define a new random variable $\varepsilon = V - U$ such that \[ f_{U,\varepsilon\vert \eta} (u,\varepsilon \vert \eta) = f_{V\vert\eta} (\varepsilon + u \vert \eta ) \left( g(Z,\delta )\right)^{-1} f_{U_0\vert \eta }( \left( g(Z,\delta )\right)^{-1} u \vert \eta). \] We can thus write
where $\tilde{\sigma}^2_U(Z) = \sigma^2_U g^2\left( Z,\delta\right) \left( 1 - \rho_{U}^\prime C^{-1}_\eta \rho_{U}\right)$, and $\tilde{\sigma}^2_V = \sigma^2_V \left( 1 - \rho_{V}^\prime C^{-1}_\eta \rho_{V} \right)$.
By tedious computations that we detail in Appendix (ref), and after integrating with respect to $U$, we obtain
with \[ \lambda(Z) = \frac{\tilde{\sigma}_U(Z)}{\tilde{\sigma}_V}, \text{ and }\sigma^2(Z) = \tilde{\sigma}^2_V + \tilde{\sigma}^2_U(Z), \] and $\Phi$ the cdf of a standard normal distribution. The distribution of $\varepsilon$ given $\eta$ is a mixture of two conditional extended skew-normal distributions, where the mixing probabilities depend on $\rho_U$ azzalini2013. When $\rho_U=\mathbf{0}$, that is, all the elements of $\rho_U$ are equal to zero, the mixing probabilities are both equal to $0.5$, and the conditional distribution of $\varepsilon$ reduces to a skew-normal distribution. That is, our specification reduces to a stochastic frontier model where the regressors are independent of the inefficiency term $U_0$.
The full information likelihood function is therefore given by \[ \mathcal{L}(\theta) = f_{\varepsilon \vert \eta} (\varepsilon \vert \eta ; \beta, \delta, \sigma^2_V,\sigma^2_U,\rho_V,\rho_U) f_{\eta}(\eta; \gamma, diag(D_\eta), ve(C_\eta) ), \] where $\theta = (\beta^\prime,\delta^\prime,\sigma^2_V,\sigma^2_U,\rho^\prime_V,\rho^\prime_U,\gamma^\prime,diag(D_\eta)^\prime, ve(C_\eta)^\prime)^\prime$; $diag(D_\eta)$ denotes the diagonal of the matrix $D_\eta$; and $ve(\cdot)$ denotes the half-vectorization of the matrix $C_\eta$ which only keeps the $(p+k)(p+k-1)/2$ elements below the main diagonal (as the matrix is symmetric and the elements of the main diagonal are equal to $1$ by construction).\footnote{This operation is defined more formally as $ve(C_\eta) = L vech(C_\eta)$, where $L$ is an elimination matrix of dimension $(p+k)(p+k-1)/2 \times (p+k)(p+k+1)/2$, which only keeps the off-diagonal elements of the half-vectorization of the matrix $C_\eta$.}
Let $\ell(\theta) = \log \mathcal{L}(\theta)$ be the log-likelihood function, and assume that $E \left[ \vert \ell(\theta) \vert \right] < \infty$ for all $\theta \in \Theta$. We define
As we can restrict $\Theta$ to be a compact parameter space, and the likelihood function is continuous in $\theta$, there exists a parameter vector $\theta_0$ which satisfies (ref) gourieroux1990.
We focus our identification analysis on the parameter $\rho_U$. To this end, we maintain the following assumption.
This assumption imposes that the parameter $\theta_1$ is first-order locally identified sargan1983. In particular, we require that the variance of the inefficiency term $\sigma^2_{U,0} > 0$. lee1986 and lee1993 have shown that when $\sigma^2_{U,0} = 0$, the stochastic frontier model is not first-order identified. Moreover, in our model, whenever $\sigma^2_U = 0$, $(\delta,\rho_U)$ are not identified. We believe this case is worthy of future investigation, but we rule it out here for simplicity.
Part (i) states that if $\rho_{U,0}$ is a solution of the maximization problem in (ref), so is $-\rho_{U,0}$, where the negative sign is applied to all components of the vector $\rho_{U,0}$. That is, the sign of all components of $\rho_U$ is not identified. Part (ii) implies that $\rho_{U,0} = \mathbf{0}$ is always a solution of (ref), which is true for any value of $\theta_1$. This result entails that the matrix of second derivatives has rank equal to $dim(\theta) - p - k$, and the model is not first-order identified at $\rho_U = \mathbf{0}$. In Part (iii), we show that first-order identification is restored when at least one component of $\rho_{U,0}$ is non-zero. A proof of this Proposition is provided in Appendix (ref).\footnote{A similar identification problem arises in Zero Inefficiency Stochastic Frontier models, see kumbhakar2013,rho2015.}
Figure (ref) illustrates the result of Part (i) of Proposition (ref). In this example, there are two endogenous regressors, one among the inputs and one among the environmental variables, so that $p = k = 1$, and the true value of $\rho_U = (0.5,0.5)^\prime$. The solid black lines are the level curves of the log-likelihood as a function of $\rho_U$, when all other parameters are taken to be known. The red dots designate the points where the log-likelihood function reaches its maximum. We can observe how both $(-0.5,-0.5)^\prime$ and $(0.5,0.5)^\prime$ are maxima of the log-likelihood function. Moreover, it can be seen from the level curves that, in this case, the log-likelihood also has a local minimum at $\rho_U = (0,0)^\prime$.
We deal with the lack of identification of the sign of $\rho_U$ by restricting the support of one of its components, say the first one, $\rho_{U1}$, to be $[0,1]$ sundberg1974. The sign of all other components is identified relative to this normalization, and provided $\rho_{U1,0}$ is in the interior of $[0,1]$.\footnote{An alternative approach would be to construct confidence sets for the identified set following chen2018, but we do not pursue it here.} When $\rho_{U1} = 0$, the sign of the other components of $\rho_U$ remains unidentified. We can therefore distinguish two cases. In the first case, there is at least one component of $\rho_U$ which is non-zero. That is, $\rho_{U1,0} > 0$, and all other components are in the interior of $[-1,1]^{p + k - 1}$. In this case, the model is first order identified, and all the components of $\rho_U$ are identified up to a sign normalization (see Proposition (ref)(iii)). In the second case, $\rho_{U} = \mathbf{0}$, and, because of Proposition (ref)(ii), the model is not first-order identified. We refer interested readers to Section (ref), where we informally discuss the choice of $\rho_{U1}$ in practice.
We let $\bar{\Theta}$ to be the parameter's space which embeds the restriction on $\rho_{U1}$, and we redefine
which exists and is (locally) unique under Assumption (ref).
While Part (ii) of Proposition (ref) states that the model is not first-order identified when $\rho_{U,0} = \mathbf{0}$, in the next proposition we show that the model is second-order identified at $\rho_{U,0} = \mathbf{0}$.
A proof is provided in Appendix (ref).
We consider an iid sample drawn from the joint distribution of $(Y,X,Z,W)$, that we denote $\lbrace (Y_i,X_i,Z_i,W_i),i = 1,\dots,n\rbrace$, where each observation follows the model in equation (ref).
Estimation is straightforward and follows from the specification of the likelihood function derived above. For all $i = 1,\dots,n$, we can write
with $\eta_i = (\eta^\prime_{X,i},\eta^\prime_{Z,i})^\prime$ and
Letting, $\ell_n(\theta) = \log \mathcal{L}_n(\theta)$ to be the sample log-likelihood function, we have \[ \hat\theta_n = \operatorname*{arg\,max}_{\theta \in \bar{\Theta}} \ell_n(\theta). \] In parallel with our identification study, we analyze our estimator's asymptotic properties depending on the true value of the parameter $\rho_U$.
Upon the additional assumption that $E\left[ \sup_{\theta \in \bar{\Theta}} \vert \ell(\theta) \vert \right] < \infty$, the log-likelihood function satisfies the required conditions for consistency newey1994h. We thus have that \[ \hat{\theta}_{n} \xrightarrow{p} \theta_{0}. \] Moreover, the log-likelihood function is at least twice continuously differentiable with respect to the parameter $\theta_{0}$. When $\rho_U$ is in the interior of $[0,1] \times [-1,1]^{p+k-1}$ and Assumption (ref) holds, standard theory of maximum likelihood estimation applies, and we can claim that \[ \sqrt{n} \left( \hat{\theta}_{n} - \theta_{0} \right) \xrightarrow{d} N\left(0,\mathcal{I}_{\theta_{0}}^{-1} \right), \] where $\mathcal{I}_{\theta_{0}}$ is the Fisher's information matrix.
However, the asymptotic distribution and the rate of convergence of our estimator are non-standard when all the components of $\rho_{U,0}$ are equal to $0$. In this case, it follows from the result of Proposition (ref) that we have a singular Hessian matrix, and one of the parameters of interest is at the boundary of the parameter space. This implies that we do not have the standard $\sqrt{n}$-rate of convergence, and that our estimator is not asymptotically normal sundberg1974b,andrews1999,rot2000. However, Proposition (ref) also implies that a reparametrization of the log-likelihood function allows us to obtain the rate of convergence and asymptotic distribution of our estimator.
Let $vec(\rho_U \rho_U^\prime)$ be the $(p+k)^2$ vectorization of the matrix $\rho_U \rho_U^\prime$. The following theorem gives the asymptotic properties of our estimator when $\rho_{U,0} = \mathbf{0}$.
The vector $\hat\tau_{\rho_U\rho^\prime_U}$ is the projection of a normal random vector onto $\mathrm{T}_1$ with respect to the Euclidean norm weighted by the matrix $\mathcal{I}_{\rho_U\rho^\prime_U} - \mathcal{I}_{\rho_U\rho^\prime_U\theta_1} \mathcal{I}^{-1}_{\theta_1} \mathcal{I}_{\theta_1\rho_U\rho^\prime_U}$ chernoff1954,andrews1999,rot2000. When $\rho_U$ is a scalar, $\hat\tau_{\rho_U\rho^\prime_U} = \max \lbrace Z_{\rho_U \rho^\prime_U},0\rbrace$. However, it is more cumbersome to derive the distribution of $\hat\tau_{\rho_U\rho^\prime_U}$ in closed form when the dimension of $\rho_U$ is greater than one and when there is dependence between the components of $vec(\rho_U\rho^\prime_U)$.\footnote{We show in a Supplementary Appendix that the off-diagonal elements of the $(p+k)^2 \times (p+k)^2$ matrix of fourth derivatives wrt $\rho_U$ are not zero in general.} The result in Part (ii) is a direct consequence of Part (i). However, it is worth highlighting that our estimator has rates of convergence slower than $\sqrt{n}$, and may not be asymptotically normal, depending on the true value of the parameter $\rho_U$. These results have important implications for obtaining standard errors and conducting inference on $\rho_U$.
To obtain standard errors and confidence intervals, we advocate the use of the subsampling method of politis1994; or the $m$-out-of-$n$ bootstrap of andrews2000. The subsampling method of politis1994 is consistent whenever the estimator has some asymptotic distributions (not necessarily normal) and when the rate of convergence is slower than $\sqrt{n}$.\footnote{One potential issue with the subsampling method is that one has to know the rate of convergence of the estimator. In practice, one can test first whether $\rho_U = \mathbf{0}$, and then apply the appropriate rate of convergence. Also, bertail1999 extend the subsampling method to the case when rates of convergence are unknown. We do not explore it here.} andrews2000 has shown that the $m$-out-of-$n$ bootstrap is consistent when the true parameter is at the boundaries but only when the estimator is $\sqrt{n}$ convergent. Provided that $m^2/n = o(1)$, the $m$-out-of-$n$ bootstrap is consistent when rates of convergence are slower than $\sqrt{n}$ bertail1999. We refer interested readers to andrews2010 for a recent study of asymptotic uniformity of subsampling and of the m-out-of-n bootstrap.
Furthermore, one may wish to conduct inference on the parameter $\rho_U$. In particular, a simple hypothesis to be tested is whether $X$ and $Z$ are independent of the inefficiency term, i.e. $\rho_U =\mathbf{0}$. The trinity of tests is an obvious candidate, but the implementation of these tests is not straightforward because of the non-standard asymptotic properties of the estimator of $\rho_U$.
andrews2001 studies the properties of the trinity of test when some parameters are at the boundary, although the author does not consider the issue of singularity of the Hessian matrix. His theoretical results about the Likelihood Ratio (LR) test can nonetheless be used following Theorem (ref), provided one can obtain an estimator of the information matrix under $H_0: vec(\rho_U \rho^\prime_U)= \mathbf{0}$. The critical values from the asymptotic distribution of the LR statistic are obtained by random draws from the vector $(Z_{\theta_1},Z_{\rho_U \rho^\prime_U})$, and by solving a quadratic programming problem andrews1999,andrews2001.
One important remark is about the Score test. Irrespective of the true value of $\rho_U$, the Score test has no power around $\rho_U=\mathbf{0}$. This is because, as shown in Proposition (ref), the score is always identically zero at that point.
We leave a thorough theoretical exploration of the properties of the Trinity of tests in this model for future work, but we explore some of the finite sample properties of the LR test in simulations.
Our framework is completed by an estimator of technical efficiency, $TE = \exp(-U_i)$, which is obtained from the conditional distribution of $U$ given $\varepsilon$ and $\eta$.
Let
where we have removed the dependence of $\lambda_{\star}$, $\sigma_{\star}$, $\mu_{1\star}$ and $\mu_{2\star}$ on $Z$ for simplicity. By equation (ref), we have that the joint density of $(\varepsilon,\eta)$ can be written as
The conditional density of $U$ given $\varepsilon$ and $\eta$ is then equal to
When both $U_0$ and $V$ are independent of $\eta$, this conditional density reduces to the one derived in jondrow1982.
Hence
By the properties of the cdf of the univariate normal distribution, this expression is shown to be equal to
This formula generalizes battese1988 formula for technical efficiencies to the endogenous case. Finally, the mean technical efficiency can be obtained as \[ E \left[ \exp(-U) \right] = E \left[E \left[ \exp(-U) \vert \varepsilon, \eta \right] \right], \] by the law of iterated expectations lee1978.
We replicate the simulation scheme in amsler2017. We consider the following model \[ Y_i = \beta_0 + X_{1i} \beta_1 + X_{2i} \beta_2 + V_i - U_{0i} \exp\left( Z_{1i} \delta_1 + Z_{2i} \delta_2\right), \] with $\beta_0 =\delta_1 = \delta_2 =0$ and $\beta_1 = \beta_2 =0.661$, and where the random variables $(X_{1i},Z_{1i})$ are taken to be exogenous (i.e. independent of the composite error term), and $(X_{2i},Z_{2i})$ are instead endogenous. We consider two instruments $(W_{1i},W_{2i})$, also independent of the error term.
The exogenous variables are generated independently from a normal distribution with means equal to $0$ and variances equal to $1$. These variables are equicorrelated, with correlation parameter equal to $0.5$.
We generate the triplet $(V,\eta_X,\eta_Z)$ from the following normal distribution \[
\sim N \left(
,
\right), \] with $\rho_V = (0.5,0.5)^\prime$, \[ \Sigma_\eta = C_\eta=
=
, \] and
with $\gamma = 0.316$.
We finally generate \[ U^\ast_0 \sim N \left( \sigma_U \rho^\prime_U C^{-1}_{\eta} \eta, \sigma^2_U (1 - \rho^\prime_U C^{-1}_{\eta} \rho_U) \right), \] with the stochastic inefficiency term given by $U_0 = \vert U^\ast_0 \vert$.
We consider two simulation schemes that differ because of the value of the parameter $\rho_U$. In Setting 1, we take $U_0$ to be independent of $\eta$ amsler2017. In Setting 2, we take $\rho_U = (0.5,0.5)^\prime$. In both settings, we impose that the first component of $\rho_U$ belongs to $[0,1]$. We take increasing sample sizes $n =\lbrace 250,500,1000 \rbrace$, and run $1000$ replications for each scenario.
Our estimation procedure is based on the maximization of the full likelihood in equation (ref).
There are two main issues for practical implementation of these models. First, the parameter space is often very large. To reduce the dimensionality of the optimization problem, one can first estimate the vector of parameters $\left( \gamma,diag(D_\eta),ve(C_\eta)\right)$ by OLS. Given $\left( \gamma,diag(D_\eta),ve(C_\eta)\right)$, one can then maximize the full likelihood with respect to the other parameters.
Moreover, the starting values for the remaining parameters need to be appropriately chosen, especially in nonlinear, high dimensional optimization problems like ours. To this end, we use the method of moments. We can write \[ E\left[ Y_i \vert X_i,Z_i,\eta_i\right] = \beta_0 + X_{1i} \beta_1 + X_{2i} \beta_2 + E\left[ V_i \vert \eta_i \right] - E\left[ U_{0i} \vert \eta_i \right]\exp\left( Z_{1i} \delta_1 + Z_{2i} \delta_2\right), \] using the assumption that $(U_0,V)$ is independent of $(X_2,Z_2)$ given $\eta$, with
We report both the average standard errors obtained by evaluating numerically the Hessian matrix of the full likelihood (Av. SE), the coverage of Wald-type confidence intervals (CI), and the coverage of confidence intervals obtained by the random subsampling method of politis1994 (CI$^\ast$). The nominal size for both is equal to $95\%$. The size of each subsample, $b$, should be selected in such a way that $b \rightarrow \infty$ and $b/n = o(1)$. We use $b = \lfloor n^{0.95}/\log(n)\rfloor$, where $\lfloor \cdot \rfloor$ denotes the integer part of a number. Choosing the size of each subsample in a data-driven way is still an open question and we do not explore it here politis1999. Based on our theoretical results, we expect that the inversion of the Hessian matrix is not going to provide reliable estimates of the standard errors in Setting 1.
Tables (ref) and (ref) below contain the results of these simulations. Table (ref) should be compared with Table 4, p. 138 of amsler2017. The mean and the standard deviation for most of the parameters are comparable with theirs. However, we achieve much better precision in estimating the variance of the inefficiency term, which, as indicated by amsler2017, is estimated imprecisely using the copula method. Both the bias and the variance decrease as the sample size $n$ increases, which ought to be expected from our MLE. The average standard errors computed using the inverse of the numerical Hessian are generally larger than the sampling standard deviation. However, Wald-type confidence intervals have good coverage, with the exception of those for the dependent parameter $\rho_U$ whose coverage is well below the nominal one. Subsampling confidence intervals have good coverage for the parameters of the stochastic frontier model. However, they undercover the first stage parameters, especially the variances of the control functions.
Results in Setting 2 are comparable to the results obtained above. It is worth noticing that, in line with our theory, standard errors are now estimated more precisely using the inverse of the Hessian matrix, and the coverage of Wald-type confidence intervals is much closer to the nominal one, for all parameters of the model. Subsampling confidence intervals perform similarly as above.
Finally, we discuss some simulation evidence about the LR test in this setting. For both simulation schemes, we test the composite nulls that $\rho_U= \mathbf{0}$ and $\rho_U = (0.5,0.5)^\prime$, respectively. When testing for $\rho_U= \mathbf{0}$, the critical values are approximated by simulations, whereas for $\rho_U = (0.5,0.5)^\prime$, the critical values are obtained from a $\chi^2_2$. To obtain the critical values in the former case, we first numerically approximate of the score vector at $\hat\theta_{1,n}$ at each sample point. We then stack to it the closed-form expression of the (vectorized) Hessian matrix for $\rho_U$, which is relatively straightforward to estimate. Its expression is given in the proof of Proposition (ref). Finally, we compute the sample information matrix by taking the inner product of the augmented score matrix. Using the generalized inverse of the information matrix, we simulate $10000$ values from the distribution of $(Z_{\theta_1},Z_{\rho_U\rho_U^\prime})$, and we obtain an estimator of $\hat{\tau}_{\rho_U\rho_U^\prime}$ by a weighted projection of the realizations of $Z_{\rho_U\rho_U^\prime}$ into the positive orthant. This last step is performed through quadratic programming, as explained in andrews1999,andrews2001.
Table (ref) contains the size of the LR test. The nominal sizes are $\lbrace 10\%, 5\%,1\%\rbrace$, respectively. The columns indicate the true value of $\rho_U$ used in the simulation exercise and the null hypothesis of the test. For $\rho_U = \mathbf{0}$, the test has size close to the nominal one, although it tends to be slightly conservative. When $\rho_U=(0.5,0.5)^\prime$, the LR test tends to have the opposite behavior as sizes are slightly larger than the nominal ones.
In Table (ref), we instead report the power properties of the LR test. The columns indicate the true value of $\rho_U$ used in the simulation exercise and the null hypothesis of the test. In general, the test has good power, and the power improves substantially as the sample size increases.
Finally, we report summary statistics for our estimators of technical efficiencies using the Battese-Coelli formula provided in equation (ref). To give a reference point to the reader, in both simulation schemes the marginal distribution of $U$ is a half-normal distribution with scale parameter equal to $\sigma^2_U = 2.7519$. Therefore, the true mean technical efficiency is equal to \[ E\left[ \exp(-U) \right] = 2 \exp\left( \frac{\sigma^2_U}{2}\right) \Phi \left( -\sigma_U\right) = 0.3846. \]
Our estimator gives a plausible interval for the values of technical efficiencies. The mean technical efficiency also approaches its true value as the sample size increases.
In this section, we consider an application using data on the agricultural sector in Nepal. The dataset consists of a cross-section of $600$ vegetable-cultivating farmers for the crop year 2015. The database is sourced from the International Food Policy Research Institute and the Seed Entrepreneurs' Association of Nepal (DataNepal). For more detail on the data, see IFPRI2017. The Output variable is total vegetable production measured in rupees. Land is measured as the total area cultivated in square feet. Labor is the sum of hours worked by hired laborers and the hours worked by household members. Fertilizers are the sum of organic and inorganic fertilizers, both measured in kilograms. Seeds are measured as the sum of hybrid and pollinated seeds in grams. As environmental variables we consider Education, as the proportion of household members with higher education or professional degree; \emph{Experience}, which is the number of years the farmer has been growing vegetables; and an indicator of risk diversification, \emph{Risk Div}, which is constructed as an Ogive index of relative economic diversification wasylenko1978. It is defined as
where $s_{i}$ is the proportion of land devoted by the farmer to crop $i$, $\bar{s}_i$ is the average sample proportion of land devoted to crop $i$, and $N_C$ is the total number of crops cultivated by each farmer. A higher value of the Risk Div index implies lower risk diversification. After removing missing values, we obtain a final sample of $497$ observations. Summary statistics of the variables used in the analysis are provided in Appendix (ref).
The model we estimate is the following \[ Y = X\beta + V - U_0 \exp(Z\delta), \] where
We allow for endogeneity of three inputs (Labor, Fertilizers, and Seeds) and one environmental variable (Risk Div). As instruments, we use a dummy for whether the farmer has suffered any natural shocks in the two years prior to the survey (Natural Shocks); the average years of experience of nearby farmers, as a measure of spillover effects (Peers Experience), and its square; three variables measuring the proportion of seeds that are owned by the farmer (\emph{Own Supplier}), obtained through formal channels such as an input retailer, a private seed company or representative, a government extension service or a research institute (\emph{Formal Supplier}), or informal channels such as a family member, a farmer's cooperative, gifted from a nearby farmer, friend or farmer from other villages, or landlord (\emph{Informal Supplier}); and interaction terms between these variables. \footnote{We have checked for weak instruments using the Cragg–Donald statistic, $CG_n$ cragg1993,stock2005. We obtain a value of $CG_n = 4.891$. Our instruments appear to be sufficiently strong based on the critical values reported in Table 5.4 of stock2005, with $2$ endogenous variables, $10$ instruments and a maximum size distortion between $10\%$ and $15\%$. This conclusion is speculative, as we do not know what the distribution of our estimator is under weak-instrument asymptotic. We have also computed the value of the Cragg–Donald statistic in our simulation study with two instruments and two endogenous variables and $N = 500$. The $95$ percentile of the test statistic is equal to $1.919$, which seems to support our conclusion.}
In this application, we normalize the dependence parameter between Fertilizers and $U_0$, $\rho_{U,\eta_{Fertilizer}}$, to be positive. The motivation for this choice is that the use of Fertilizers for production may be related to soil quality, unobserved by the econometrician and which ultimately influences the efficiency of the producer. We have run some robustness checks and our results are not sensitive to the choice of normalization.
Table (ref) reports the estimated coefficients and the $95\%$ confidence intervals for our empirical example. The confidence intervals are obtained using $496$ subsamples of size $b = 50$.
The left panel shows the estimation results assuming exogeneity. All the estimated coefficients for inputs are positive and significant, although generally small in magnitude, being Labor and Land, the inputs with the largest estimated effect. None of the environmental variables appears to have a significant effect on inefficiency. The exogenous model detects very little inefficiency, and the parameter $\delta$ is estimated imprecisely, as it can be noticed by the length of the confidence intervals.
In the center panel, we report the estimation results controlling for endogeneity but restricting $\rho_U = \mathbf{0}$, i.e., imposing independence between the endogenous variables and the inefficiency term. As it is usually the case in instrumental variable models, confidence intervals are wider than in the model assuming exogeneity. However, controlling for endogeneity substantially changes the conclusions obtained from this empirical example, regarding the effect of inputs, but especially the effect of the inefficiency determinants. We find that most of the estimated coefficients for the inputs are positive but significant only for Fertilizers and Seeds. Regarding the environmental variables, the estimated coefficient of Risk Div is negative, which means that farmers cultivating fewer crops (i.e., with lower risk diversification) are more efficient, possibly due to higher specialization levels. However, the effect is not significant at the $5\%$ level. The correlation between the two-sided error term and the endogenous variables is negative and significant, except for Risk Div.
The right panel shows the results controlling for endogeneity without restricting the dependence between the endogenous variables and the inefficiency term. The estimated coefficients are quite similar to those when we impose that $\rho_U = 0$, except for Experience whose coefficient is positive and significant in the third model. This result may be due to more experienced, i.e. older, farmers being more conservative, and therefore, less willing to implement new practices that could improve their efficiency coelli1996. The $95\%$ confidence intervals are much narrower than those for the restricted estimator and indicate that only the choice of Fertilizers may be related to the producer's inefficiency level. We test for the absence of dependence between the endogenous variables and the inefficiency term being equal to $0$ using the LR test, as explained in the simulation study. The value of the test statistic is equal to $6.77$, and we obtain a $95\%$ critical value equal to $11.41$. Hence, we cannot reject the null that $\rho_U = \mathbf{0}$.
We also test the null that $\sigma^2_U=0$ in the three models. Under the null, $\rho_U$ and $\delta$ are nuisance parameters. As we do not know the asymptotic distribution of the LR test statistics in this case and we are testing a parameter at the boundary, critical values are obtained from an equal mixture of a mass point at $0$ and a $\chi^2$-distribution with $1$ degree of freedom lee1993,andrews2001,ketz2018. In both the exogenous and the endogenous model, we reject the null of no inefficiency at the $1\%$ level. The value of the LR statistic is $4.66$ for the exogenous model, and $79.46$ for the endogenous model, with a critical value equal to $3.82$.
Finally, Figure (ref) reports the kernel density estimator of technical efficiency for the endogenous models. The density of the efficiency scores is similar in both models, which is an expected result given that we cannot reject the null that $\rho_U$ is equal to 0.
We propose and study an estimator of stochastic frontier models when both the production inputs and the environmental variables are correlated with the two-sided stochastic error term and the one-sided stochastic inefficiency term. Our identification and estimation strategy is based on control functions that fully capture the dependence between regressors and unobservables. While the joint density of the two-sided stochastic error term and the control function is modeled as a normal distribution, one of the main challenges for direct maximum likelihood estimation is to write the joint density of the stochastic inefficiency term and the control function in closed form. To circumvent this issue, amsler2017 use copula functions to model the dependence between observables and unobservables components of the model, and employ a simulated maximum likelihood procedure to obtain the parameter's estimate. This estimator may not be easy to implement and may be computationally slow. Moreover, instrumental variable methods lead to lower precision in the estimate and simulated methods can increase this lack of precision even further.
In this work, we provide a simple maximum likelihood estimator that aims at avoiding these potential pitfalls. Our main assumption is that the conditional distribution of the stochastic inefficiency term given the control functions is a folded normal distribution. This distribution reduces to the half-normal when there is no endogeneity. This makes our model a straightforward extension of the normal-half-normal model to include endogenous regressors. We shed light on some new identification issues, and we provide some theoretical results on estimation and inference. Our estimator is easy and fast to implement, and enjoys good finite sample properties.
Additional research on the properties of the trinity of tests and on testing the distributional assumptions on the error term is needed. Moreover, extensions of our model to panel data with time-varying endogeneity and true fixed effects could be of interest.