---
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
