spefa: spatial stochastic frontier with fixed effects and endogeneity

1. What this package estimates

spefa implements the maximum-likelihood estimator developed in

Giannini, M. (2025). A spatial stochastic frontier model with fixed effects and endogenous environmental variables. Spatial Economic Analysis, 20(3), 420–441. doi:10.1080/17421772.2024.2414962

The model is a stochastic production (or cost) frontier on panel data that simultaneously accounts for three features that the earlier literature treated separately:

• individual fixed effects \(\alpha_i\), removed by first differencing (following Wang and Ho, 2010), which avoids the incidental-parameters problem and sidesteps the singularity of the within transformation in spatial panels;
• endogenous regressors and environmental variables, corrected through a Gaussian control function obtained from a Cholesky factorisation of the error covariance (following Kutlu, Tran and Tsionas, 2020);
• spatial spillovers, through a spatial-autoregressive (SAR) term \(\lambda W y\).

For unit \(i\) and period \(t\) the frontier is

\[ y_{it} = \alpha_i + \lambda \sum_j w_{ij} y_{jt} + x_{it}'\beta - u_{it} + v_{it}, \qquad u_{it} = h_{it}\,u_i^{*}, \]

with \(h_{it} = \exp(z_{it}'\phi) > 0\) a multiplicative scaling and a time-invariant \(u_i^{*} \sim N^{+}(\mu,\sigma_u^2)\) (half-normal when \(\mu = 0\)). Endogenous variables \(q_j\) obey a reduced form \(q_{j,it} = z_{j,it}'\delta_j + \varepsilon_{j,it}\), and endogeneity is the correlation \(\rho_j = \mathrm{corr}(\varepsilon_j, v)\). The full derivation of the likelihood is in Giannini (2025); the package implements it in closed form and maximises it numerically.

The implementation depends only on base R (the stats package).

2. Installation

install.packages("spefa_0.1.0.tar.gz", repos = NULL, type = "source")
library(spefa)

3. The data you supply

The estimator needs a balanced panel and a spatial weight matrix:

• a data.frame in long format with a unit identifier and a time identifier;
• the outcome y and the frontier regressors X (these may include endogenous variables);
• for each endogenous variable, one or more instruments;
• the variables that enter the inefficiency scaling \(h_{it}\);
• an \(N \times N\) weight matrix W, with rows and columns ordered by the sorted unit id, and (typically) row-normalised.

4. The spefa() function and its options

spefa(formula, data, index, W,
      endogenous = NULL, scaling = NULL, mu = FALSE,
      control = list(maxit = 500, reltol = 1e-9))

5. What is returned, and the available methods

spefa() returns an object of class "spefa". The estimated parameter vector contains the frontier slopes \(\beta\), the instrument slopes \(\delta_j\), the scaling coefficients \(\phi\), the reduced-form variances \(\sigma^2_{\varepsilon_j}\), the endogeneity correlations \(\rho_j\), \(\sigma^2_v\), \(\sigma^2_u\), (optionally \(\mu\)) and the spatial parameter \(\lambda\).

Standard errors come from the Hessian at the optimum (computed on the unconstrained scale) propagated to the natural parameters by the delta method.

6. A worked example

A small demo panel (spefademo, \(N = 40\), \(T = 12\)) ships with the package and is used here. The data-generating values are \(\beta=(0.5,0.5)\), \(\delta=(1,1)\), \(\phi=(0.3,0.3)\), \(\sigma^2_\varepsilon=0.09\), \(\sigma^2_v=0.04\), \(\rho=(0.5,0.5)\), \(\sigma^2_u=0.09\), \(\lambda=0.5\); q1 and q2 are endogenous, q1 enters the frontier and q2 the scaling.

library(spefa)
data(spefademo)

fit <- spefa(y ~ x1 + q1, data = spefademo, index = c("id", "time"),
             W = spefaW, endogenous = list(q1 = ~ z1, q2 = ~ z2),
             scaling = ~ x2 + q2)

summary(fit)

Indicative output from a larger sample (\(N=150\), to show parameter recovery; the bundled \(N=40\) panel gives the same point estimates with wider standard errors):

               Estimate Std.Error  z value  Pr(>|z|)
x1              0.5002   0.0027    182.2    <2e-16 ***
q1              0.5007   0.0032    155.2    <2e-16 ***
delta.q1.z1     0.9991   0.0056    179.8    <2e-16 ***
delta.q2.z2     1.0029   0.0054    186.6    <2e-16 ***
x2              0.2995   0.0120     24.9    <2e-16 ***
q2              0.3051   0.0120     25.5    <2e-16 ***
rho.q1          0.4763   0.0134     35.5    <2e-16 ***
rho.q2          0.5095   0.0130     39.2    <2e-16 ***
sigma2_v        0.0391   0.0010     38.5    <2e-16 ***
sigma2_u        0.0879   0.0162      5.4    5e-08  ***
lambda          0.5027   0.0044    114.5    <2e-16 ***

Spatial parameter lambda = 0.5027
Endogeneity correlations rho (q1, q2): 0.476, 0.509

impacts(fit)            # direct / indirect / total marginal effects
ef <- efficiency(fit)   # spatially-corrected technical efficiency
summary(as.vector(ef$eff))

7. Reading the results

• Endogeneity. A \(\rho_j\) statistically different from zero signals endogeneity of variable \(j\); the \(z\)-test in the coefficient table is a direct test of \(H_0:\rho_j=0\). When all \(\rho_j=0\) the control-function correction vanishes and the estimator collapses to the exogenous case.
• Spatial dependence. \(\lambda\) measures the strength of spillovers; the impacts() table splits each regressor’s effect into a direct (own-unit) and an indirect (spillover) component, with total \(\approx \beta_k/(1-\lambda)\) under row-normalised W.
• Efficiency. efficiency(fit)$eff gives unit-by-period technical efficiency in \((0,1]\); with spatial = TRUE the scores incorporate neighbours’ inefficiency through \((I-\lambda W)^{-1}\).

8. Practical notes and identification

• The panel must be balanced, and the rows/columns of W must follow the sorted unit id used in data.
• The inefficiency parameters \(\phi\) and \(\sigma_u\) are identified only when the scaling variables vary over time (otherwise first differencing removes the inefficiency) and when the inefficiency is non-negligible. If the inefficiency is very small, expect \(\phi\) and \(\sigma_u\) to be weakly identified.
• The endogeneity correction maintains the usual control-function assumptions: joint normality of the structural and reduced-form errors, and inefficiency independent of the regressors.
• The likelihood evaluates the reduced-form (instrument) density on the data and counts the Gaussian normalising constant once per unit; the truncated-normal terms use the log-CDF for numerical stability.

References

Giannini, M. (2025). A spatial stochastic frontier model with fixed effects and endogenous environmental variables. Spatial Economic Analysis, 20(3), 420–441.

Jondrow, J., Lovell, C. A. K., Materov, I. S., & Schmidt, P. (1982). On the estimation of technical inefficiency in the stochastic frontier production function model. Journal of Econometrics, 19(2–3), 233–238.

Kutlu, L., Tran, K. C., & Tsionas, M. G. (2020). A spatial stochastic frontier model with endogenous frontier and environmental variables. European Journal of Operational Research, 286(1), 389–399.

LeSage, J., & Pace, R. K. (2009). Introduction to Spatial Econometrics. CRC Press.

Wang, H.-J., & Ho, C.-W. (2010). Estimating fixed-effect panel stochastic frontier models by model transformation. Journal of Econometrics, 157(2), 286–296.