---
title: "Replicating published results"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Replicating published results}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
```

```{r setup}
library(pgt)
```

The estimators with a numerical example printed in their source paper
are validated against it: the weak-G-disposability and directional
models against Rødseth (2025) Tables 2 and 3, and the by-production
model against the analytic Example 1 of Murty, Russell and Levkoff
(2012). These replications run as unit tests on every check; this
vignette reproduces them in the open. The remaining estimators carry
analytic hand-computed checks in the test suite rather than
published-table replications; the final section shows the
internal-consistency check for the materials-balance cost model.

## Rødseth (2025), Table 2: minimal controlled emissions

The `pigfarms` data is Table 1 of Rødseth (2025): five pig-finishing
farms with feed and piglets as material inputs, meat as the good output,
and controlled, uncontrolled and abatement emission columns for
nitrogen. Because the paper does not print the nitrogen coefficients of
feed and piglets, the uncontrolled-emission aggregate replaces them as
the single material carrier ($u = 1$, $v = 0$); the uncontrolled
emission equals controlled plus abatement, so the accounts close
exactly. Keeping all five columns as inputs with zero material
coefficients on feed and piglets (the mapping in the `pigfarms` help
page) yields the same Table 2 minima; the test suite asserts the two
mappings agree.

```{r}
data(pigfarms)
pigfarms
```

```{r}
tech <- pgt_tech(
  x = pigfarms[, c("uncontrolled", "labor", "capital")],
  y = pigfarms$meat,
  b = pigfarms$controlled,
  u = c(1, 0, 0),
  a = pigfarms$abatement,
  id = pigfarms$farm
)
fit <- pgt(tech, model = "wgd")
data.frame(farm = fit$results$id, b_star = fit$results$b_star,
           z_star = fit$results$z_star, a_star = fit$results$a_star)
```

The minimal controlled emissions are $(16, 16, 16, 20, 16)$, which match
the paper's Table 2. Because the abatement output is observed, the fit
also reports the projection's implied uncontrolled emission `z_star`
and abatement `a_star`; for farms A, D and E these match the paper's
Table 2 uncontrolled and abatement entries, while for the frontier
farms the optimal peer mix is not unique and the split between
`z_star` and `a_star` is solver-dependent.

## Rødseth (2025), Table 3: the directional representation

The factorially determined multi-output model (Eq. 13, abatement fixed at
each farm's own level) jointly expands the good output and contracts the
bad output. Its good-output coefficient $v = 7/6$ follows from farms A
and B, which share inputs: their meat differs by 3 and their uncontrolled
emissions by 3.5. Each farm's feed coefficient closes its account to the
uncontrolled column.

```{r}
v <- 7 / 6
u_feed <- (pigfarms$uncontrolled + v * pigfarms$meat) / pigfarms$feed
tech3 <- pgt_tech(
  x = pigfarms[, c("feed", "piglet", "labor", "capital")],
  y = pigfarms$meat,
  b = pigfarms$controlled,
  u = cbind(u_feed, 0, 0, 0),
  v = v,
  a = pigfarms$abatement,
  id = pigfarms$farm
)
fdmo <- pgt(tech3, model = "fdmo")
r3 <- fdmo$results
comparison <- data.frame(
  farm = r3$id,
  good_eff = round(r3$good_eff, 3),
  paper_good = c(0, 0, 3, 0, 0.8),
  bad_eff = round(r3$bad_eff, 3),
  paper_bad = c(0, 0, 3.5, 0, 1.0),
  maximal_y = round(r3$maximal_y, 3)
)
comparison
```

Farms A, B and D are on the frontier (gross inefficiency 0), and farm C
expands its good output by 3 and cuts its bad output by 3.5, a gross
inefficiency of 6.5: farms A to D reproduce Table 3 exactly. Farm E is
the one divergence, shown side by side above. The paper prints
$(\theta_y, \theta_b) = (0.8, 1.0)$; the package returns
$(0.789, 0.921)$. The two cannot both be exact: whenever the
uncontrolled emission equals controlled plus abatement, the
materials-balance row forces $\theta_b = v\,\theta_y$ for every farm,
and $7/6 \times 0.8 \approx 0.93$, not $1.0$. The package's pair
satisfies the forced ratio to machine precision, so we read the printed
row as rounded; farm E's maximal good output still rounds to the
paper's 10.8.

## Murty, Russell and Levkoff (2012): the by-production example

Their Example 1 is a constant-returns technology with one input, one
good output $y$ and one bad output, written $z$ in the paper (our $b$),
solved analytically in the paper.

```{r}
mrl <- pgt_tech(
  x = matrix(c(1, 1, 1, 2, 2), ncol = 1, dimnames = list(NULL, "x")),
  y = c(2, 3 / 2, 2 / 3, 3, 2),
  b = c(4, 1, 2, 5, 3),
  u = 1, polluting = "x",
  id = paste0("DMU", 1:5)
)
bp <- pgt(mrl, model = "byprod", returns = "crs")
rb <- bp$results
data.frame(id = rb$id, output_eff = round(rb$output_eff, 4),
           emission_eff = round(rb$efficiency, 4), mean_eff = round(rb$mean_eff, 4))
```

For DMU3 the output efficiency is $1/3$, the emission efficiency $1/2$
and the graph measure $5/12$; for DMU2 they are $3/4$, $1$ and $7/8$.
These match the analytic scores in the paper.

## Coelli, Lauwers and Van Huylenbroeck (2007): coefficients and the EE = TE x EAE decomposition

This section is an internal-consistency check, not a replication of a
printed table: the paper's farm-level data are not published, so the
four farms below are invented. What comes from the paper are the
phosphorus contents (their footnote 16): feed 0.0124, piglet 0.0117 and
meat 0.0117 kg P per kg. The check confirms that the implemented cost
model decomposes environmental efficiency as $EE = TE \times EAE$ with
these coefficients; the model's LPs are verified against hand-computed
solutions in the test suite.

```{r}
farms <- data.frame(
  feed = c(200, 210, 190, 205),
  piglet = c(20, 22, 18, 21),
  meat = c(100, 100, 95, 102)
)
tech_c <- pgt_tech(
  x = as.matrix(farms[, c("feed", "piglet")]),
  y = farms$meat,
  b = 0.0124 * farms$feed + 0.0117 * farms$piglet - 0.0117 * farms$meat,
  u = c(feed = 0.0124, piglet = 0.0117),
  v = 0.0117
)
mbc <- pgt(tech_c, model = "mb_cost", returns = "crs")
rc <- mbc$results
data.frame(id = rc$id, EE = round(rc$efficiency, 4),
           TE = round(rc$te, 4), EAE = round(rc$eae, 4))
all.equal(rc$efficiency, rc$te * rc$eae)
```

The identity $EE = TE \times EAE$ holds to machine precision. (The
identity holds by construction of the reported components, so this
check guards the wiring, not the mathematics; the LP solutions
themselves are checked by hand in the test suite.)

## References

- Coelli, T., Lauwers, L., & Van Huylenbroeck, G. (2007). Environmental
  efficiency measurement and the materials balance condition. *Journal
  of Productivity Analysis*, 28(1-2), 3-12. doi:10.1007/s11123-007-0052-8
- Murty, S., Russell, R. R., & Levkoff, S. B. (2012). On modeling
  pollution-generating technologies. *Journal of Environmental
  Economics and Management*, 64(1), 117-135.
  doi:10.1016/j.jeem.2012.02.005
- Rødseth, K. L. (2025). On the development of a unified, nonparametric
  materials balance-based efficiency analysis model and its
  applications. *Journal of Productivity Analysis*, 64(3), 305-319.
  doi:10.1007/s11123-025-00768-0
