--- title: "Productivity change and inference" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Productivity change and inference} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 6, fig.height = 4) ``` ```{r setup} library(pgt) ``` ## Productivity change with a bad output 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 output distance function, the output counterpart of the directional input distance function of 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; the global index of Oh (2010) pools every period into one reference technology, which removes the infeasibility and makes the index circular. It does not address the separate inconsistency shown by Aparicio, Pastor and Zofío (2013), in which the technical-change component signals regress although the frontier moves toward more good and less bad output. `pgt_ml()` computes the global index under the pollution-generating technology and splits it into efficiency change and best-practice gap 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. The default `"wd"` imposes weak disposability of the bad output in the Kuosmanen (2005) form. With `returns = "crs"` its abatement weights are redundant and, for one bad output, it is the constant-returns technology of Chung, Färe and Grosskopf (1997, Eq. 3.9), which Oh (2010) also uses; both papers index several bad outputs jointly, whereas `pgt_ml()` indexes one pollutant at a time; the variable-returns default extends it. The `"input_envelope"` and `"envelope"` options treat the bad output as reducible to the peer emission envelope and are exploratory companions; the input-keeping variant matches the shape of the input-fixed benchmark of `pgt()` rather than its Eq. 6 programme. ### A constructed panel The panel below has a frontier that improves between the two periods: the same output is produced with 15 percent less emission. ```{r} 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) ``` The median index exceeds 1, and best-practice gap change carries it: the frontier moved out while each unit kept its position relative to the frontier, so efficiency change stays near 1. ### A real rice panel The Philippine rice panel in the `frontier` package spans eight years. Treating fertiliser nitrogen as the pollutant potential gives a nitrogen-based productivity index. ```{r, eval = requireNamespace("frontier", quietly = TRUE)} 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) ``` ## Subsampling inference Efficiency 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). Subsampling gives valid intervals under minimal assumptions (Politis and Romano 1994; Politis, Romano and Wolf 1999); for DEA, Kneip, Simar and Wilson (2008) prove the m-out-of-n bootstrap, which draws with replacement, consistent under variable returns for the score of a fixed point, and Simar and Wilson (2011) note that the resampling can also be done without replacement, although in their simulations resampling with replacement gave coverage closer to nominal. `boot_pgt()` covers cross-sectional efficiency scores only; no inference is currently available for the `pgt_ml()` index. The example takes the sixty plants of one year of the synthetic `steeldemo` panel, an independent cross-section. ```{r} data(steeldemo) steel60 <- steeldemo[steeldemo$year == 2021, ] tech <- pgt_tech( x = steel60[, c("coal_coke", "other_fuel", "raw_material", "flux", "capture_energy")], y = steel60$production, b = steel60$emissions, a = steel60$captured, v = 0.01467, x_abate = "capture_energy", group = steel60$route, id = steel60$plant ) bt <- boot_pgt(tech, model = "wgd", B = 40, seed = 1) head(bt$per_dmu) ``` Each 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). ```{r} boot_pgt_sensitivity(tech, model = "wgd", B = 25, m_grid = c(15, 25, 40), seed = 1) ``` A 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). ### Coverage under a known data-generating process 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 with exponential inefficiency and two known frontiers, the affine $b_{\min}(y) = 2 + 0.5y$ and the curved $b_{\min}(y) = 1 + 0.05y^2$, where the true score $\theta = b_{\min}(y)/b$ is the estimand. Two hundred replicates per cell, `B = 100`, nominal level 95 per cent. A unit gets no interval when fewer than `B/2` of its subsample programmes are feasible, so the study reports availability, coverage among units with an interval ("given interval") and the share of all units whose interval exists and covers ("delivered"). The table shows the cells at the defaults $m = \mathrm{round}(L^{0.7})$ and $\kappa = 2/3$ from the precomputed results shipped with the package. ```{r, echo = FALSE} cov <- read.csv(system.file("simulations", "coverage-results.csv", package = "pgt")) dflt <- abs(cov$kappa - 2 / 3) < 1e-9 & cov$m == round(cov$L^0.7) knitr::kable( cov[dflt, c("dgp", "L", "available", "coverage_cond", "coverage_delivered", "coverage_frontier", "coverage_interior", "coverage_group_means")], digits = 3, row.names = FALSE, col.names = c("frontier", "L", "available", "given interval", "delivered", "frontier (theta > 0.9)", "interior", "group means") ) ``` On the affine frontier the intervals over-cover, and on the curved one they under-cover by two to three points at both sample sizes. The missing intervals fall on the units with the largest output. The shipped results also vary `m` and `kappa` at `L = 100`: a larger subsample or a faster rate narrows the intervals and lowers coverage, to 0.77 under `kappa = 1` on the curved frontier. One simulation design is evidence, not proof: the DEA rate remains borrowed for the materials-balance programmes, which is why the intervals keep their heuristic label and why `boot_pgt_sensitivity()` belongs in the workflow. ## References 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