--- 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 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. ### 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 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), 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. ```{r} 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) 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 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. ```{r, echo = FALSE} cov <- read.csv(system.file("simulations", "coverage-results.csv", package = "pgt")) knitr::kable( cov[, c("L", "reps", "coverage_all", "coverage_frontier", "coverage_interior", "coverage_group_means")], digits = 3, col.names = c("L", "replicates", "all DMUs", "frontier (theta > 0.9)", "interior", "group means") ) ``` 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. ## 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