The Malmquist-Luenberger index (Chung, Färe and Grosskopf 1997)
measures productivity change when part of the output is undesirable, by
scoring each unit on a directional distance function (Chambers, Chung
and Färe 1996) that expands the good output and contracts the bad output
together. The adjacent-period version can be infeasible across periods
(Aparicio, Pastor and Zofío 2013); the global index of Oh (2010) pools
every period into one reference technology, which removes the
infeasibility and makes the index circular. pgt_ml()
computes it under the pollution-generating technology and splits it into
efficiency change and best-practice change, \(GML = EC \times BPC\), with \(GML > 1\) signalling growth.
This estimator is experimental: the global Malmquist-Luenberger index under a materials-balance technology is not yet settled in the literature, and the default technologies do not re-impose the per-DMU materials-balance cap on the projected point.
The technology argument selects the reference
technology. technology = "wd" imposes weak disposability of
the bad output in the Kuosmanen (2005) form; this is the technology
under which Chung, Färe and Grosskopf (1997) and Oh (2010) define the
index, so it is the faithful comparator. The default "wgd"
and the "envelope" option treat the bad output as reducible
to the peer emission envelope and are exploratory companions.
The panel below has a frontier that improves between the two periods: the same output is produced with 15 percent less emission.
set.seed(1)
n <- 15
d1 <- data.frame(id = 1:n, y = runif(n, 5, 10), x = runif(n, 8, 12),
b = runif(n, 2, 6))
d2 <- data.frame(id = 1:n, y = d1$y * 1.03, x = d1$x, b = d1$b * 0.85)
d <- rbind(cbind(d1, period = 1), cbind(d2, period = 2))
panel <- pgt_tech(x = d[, "x", drop = FALSE], y = d$y, b = d$b,
period = d$period, id = d$id)
ml <- pgt_ml(panel)
summary(ml)
#> pgt global Malmquist-Luenberger index (technology = wgd, returns = vrs)
#>
#> Distribution (GML > 1 is productivity growth):
#> 0% 25% 50% 75% 100%
#> GML 1.03 1.033 1.0347 1.0499 1.1249
#> EC 1.00 1.000 1.0000 1.0000 1.0000
#> BPC 1.03 1.033 1.0347 1.0499 1.1249
#>
#> Geometric mean by period transition:
#> to n gml ec bpc
#> 2 15 1.049 1 1.049The median index exceeds 1, and best-practice change carries it: the frontier moved out while each unit kept its position relative to the frontier, so efficiency change stays near 1.
The Philippine rice panel in the frontier package spans
eight years. Treating fertiliser nitrogen as the pollutant potential
gives a nitrogen-based productivity index.
data("riceProdPhil", package = "frontier")
r <- riceProdPhil
uN <- 0.46
rice_panel <- pgt_tech(
x = as.matrix(r[, c("AREA", "LABOR", "NPK", "OTHER")]),
y = r$PROD,
b = uN * r$NPK,
u = c(AREA = 0, LABOR = 0, NPK = uN, OTHER = 0),
v = 0,
period = r$YEARDUM,
id = as.character(r$FMERCODE)
)
rice_ml <- pgt_ml(rice_panel)
summary(rice_ml)
#> pgt global Malmquist-Luenberger index (technology = wgd, returns = vrs)
#>
#> Distribution (GML > 1 is productivity growth):
#> 0% 25% 50% 75% 100%
#> GML 0.5981 0.9082 1.0000 1.1291 1.7404
#> EC 0.6339 0.9410 1.0000 1.0802 1.5929
#> BPC 0.6907 0.9203 0.9997 1.0971 1.4926
#>
#> Geometric mean by period transition:
#> to n gml ec bpc
#> 2 43 0.9809 1.0375 0.9455
#> 3 43 1.0464 1.0605 0.9866
#> 4 43 1.0046 0.9807 1.0244
#> 5 43 0.9322 0.9687 0.9623
#> 6 43 1.0069 0.9755 1.0322
#> 7 43 0.9027 0.9925 0.9095
#> 8 43 1.2128 1.0071 1.2043Efficiency scores are point estimates from a finite sample of DMUs.
boot_pgt() gauges their sampling variability by re-solving
each unit against random subsamples of the reference set (m out of n,
without replacement), the resampling scheme with the firmest footing for
boundary estimators (Politis and Romano 1994; Politis, Romano and Wolf
1999; Kneip, Simar and Wilson 2008). boot_pgt() covers
cross-sectional efficiency scores only; no inference is currently
available for the pgt_ml() index. The example takes the
first 60 plants of the synthetic steeldemo panel.
data(steeldemo)
steel60 <- steeldemo[1:60, ]
tech <- pgt_tech(
x = steel60[, c("coal_coke", "other_fuel", "raw_material", "flux")],
y = steel60$production, b = steel60$emissions, v = 0.01467,
group = steel60$route, id = steel60$plant
)
bt <- boot_pgt(tech, model = "wgd", B = 40, seed = 1)
#> Warning: 17 of 60 DMUs have infeasible LPs in more than half of their subsample
#> replicates; their intervals rest on few replicates (see per_dmu$n_ok).
head(bt$per_dmu)
#> id group estimate lower upper se n_ok
#> 1 P001 Integrated 0.7202900 0.6521703 0.7202900 0.01888103 40
#> 2 P001 Integrated 0.7171565 0.6585343 0.7171565 0.01620572 40
#> 3 P001 Integrated 0.7081322 0.6434389 0.7081322 0.01730729 40
#> 4 P002 Integrated 0.9461050 0.9078969 0.9425286 0.01503927 19
#> 5 P002 Integrated 0.8494481 0.7783673 0.8494481 0.02262132 34
#> 6 P002 Integrated 0.8856638 0.8046069 0.8856638 0.02859469 27Each unit is evaluated against the subsample frontier alone; the unit
itself is not appended to the subsample. Replicates in which its
programme is infeasible are dropped and counted in
per_dmu$n_ok, and the point estimate can lie above its own
upper endpoint, because the interval targets the true score, which the
estimate overestimates.
The intervals follow the subsampling recipe: the subsample scores are recentred at the point estimate and rescaled by the relative convergence rate \((m/L)^\kappa\), so the interval extends downward from the estimate, matching the direction of the frontier bias (a smaller reference sample can only move the estimated frontier inward). The intervals remain heuristic: resampling theory for data envelopment analysis covers the radial and directional technologies, not the materials-balance programs here, so \(\kappa\) borrows the DEA rate as an approximation. Before relying on the intervals, check that their width is stable in the subsample size (Simar and Wilson 2011).
boot_pgt_sensitivity(tech, model = "wgd", B = 25,
m_grid = c(15, 25, 40), seed = 1)
#> Warning: 22 of 60 DMUs have infeasible LPs in more than half of their subsample
#> replicates; their intervals rest on few replicates (see per_dmu$n_ok).
#> Warning: 11 of 60 DMUs have infeasible LPs in more than half of their subsample
#> replicates; their intervals rest on few replicates (see per_dmu$n_ok).
#> Warning: 1 of 60 DMUs have infeasible LPs in more than half of their subsample
#> replicates; their intervals rest on few replicates (see per_dmu$n_ok).
#> m median_width median_width_pos mean_se
#> 1 15 0.07378358 0.07655731 0.02642923
#> 2 25 0.06319221 0.06398525 0.02588980
#> 3 40 0.04998538 0.05013635 0.01991171A width that changes little across the grid supports the chosen subsample size; a strong trend means the size, not the data, is driving the interval.
When the technology carries groups, boot_pgt() also
reports intervals for the group mean scores. These are a descriptive
stability band, not supported inference: the replicate statistic
averages all group members against a common subsample frontier, so
replicate variation reflects movement of the frontier only, and the
per-DMU convergence rate is reused for the mean (Kneip, Simar and Wilson
2015 give the theory for aggregates of efficiency scores).
The heuristic label invites the obvious question: how close to
nominal do these intervals get when the truth is known? The package
ships a Monte Carlo study (inst/simulations/coverage.R) for
the envelope model under a data-generating process with a known linear
frontier \(b_{\min}(y) = 2 + 0.5y\) and
exponential inefficiency, where the true score \(\theta = b_{\min}(y)/b\) is the estimand.
Two hundred replicates per sample size, B = 100, nominal
level 95 per cent; the table reports the precomputed results shipped
with the package.
| L | replicates | all DMUs | frontier (theta > 0.9) | interior | group means |
|---|---|---|---|---|---|
| 50 | 200 | 0.943 | 0.957 | 0.937 | 0.950 |
| 100 | 200 | 0.970 | 0.980 | 0.966 | 0.978 |
Empirical coverage sits near the nominal 95 per cent for both sample sizes, for frontier and interior units alike, and errs on the conservative side at the larger sample. One simulation under one data-generating process is evidence, not proof: the DEA rate remains borrowed for the materials-balance programmes, which is why the intervals keep their heuristic label.
Aparicio, J., Pastor, J. T., & Zofío, J. L. (2013). On the inconsistency of the Malmquist-Luenberger index. European Journal of Operational Research, 229(3), 738-742. doi:10.1016/j.ejor.2013.03.031
Chambers, R. G., Chung, Y., & Färe, R. (1996). Benefit and distance functions. Journal of Economic Theory, 70(2), 407-419. doi:10.1006/jeth.1996.0096
Chung, Y. H., Färe, R., & Grosskopf, S. (1997). Productivity and undesirable outputs: A directional distance function approach. Journal of Environmental Management, 51(3), 229-240. doi:10.1006/jema.1997.0146
Kneip, A., Simar, L., & Wilson, P. W. (2008). Asymptotics and consistent bootstraps for DEA estimators in nonparametric frontier models. Econometric Theory, 24(6), 1663-1697. doi:10.1017/S0266466608080651
Kneip, A., Simar, L., & Wilson, P. W. (2015). When bias kills the variance: Central limit theorems for DEA and FDH efficiency scores. Econometric Theory, 31(2), 394-422. doi:10.1017/S0266466614000413
Kuosmanen, T. (2005). Weak disposability in nonparametric production analysis with undesirable outputs. American Journal of Agricultural Economics, 87(4), 1077-1082. doi:10.1111/j.1467-8276.2005.00788.x
Oh, D.-h. (2010). A global Malmquist-Luenberger productivity index. Journal of Productivity Analysis, 34(3), 183-197. doi:10.1007/s11123-010-0178-y
Politis, D. N., & Romano, J. P. (1994). Large sample confidence regions based on subsamples under minimal assumptions. The Annals of Statistics, 22(4), 2031-2050. doi:10.1214/aos/1176325770
Politis, D. N., Romano, J. P., & Wolf, M. (1999). Subsampling. Springer, New York. doi:10.1007/978-1-4612-1554-7
Simar, L., & Wilson, P. W. (2011). Inference by the m out of n bootstrap in nonparametric frontier models. Journal of Productivity Analysis, 36(1), 33-53. doi:10.1007/s11123-010-0200-4