Introduction to pgt

library(pgt)

The materials-balance principle

A pollution-generating technology turns material inputs into a good output and an unavoidable bad output. The pollutant content of the inputs has to end up somewhere: embodied in the product, removed by abatement, or emitted. pgt builds efficiency models on that accounting identity, \[u'x_l - v y_l \ge b_l,\] where \(x_l\) are the material inputs of unit \(l\), \(u\) their pollutant flow coefficients (for example tonnes of CO2 potential per unit of input), \(y_l\) the good output, \(v\) the pollutant retained in the product, and \(b_l\) the emissions. The gap \(u'x_l - v y_l - b_l\) is the pollutant that leaves the process other than through the measured emission. A note on terminology: the mass balance is an identity when every exit is measured; the inequality form above allows for unmeasured abatement, and when the abatement output \(a_l\) is observed the account closes as the equality \(u'x_l - v y_l = b_l + a_l\). Following Rødseth (2025), the documentation calls both forms the materials-balance identity.

The principle itself is old: Ayres and Kneese (1969) built the materials-balance view of production and externalities, Lauwers (2009) made the case for carrying it into frontier-based eco-efficiency models, and Dakpo, Jeanneaux and Latruffe (2016) survey how pollution-generating technologies have been modelled in nonparametric benchmarking since.

Existing R packages handle undesirable outputs by transforming the data (the Seiford and Zhu 2002 translation in deaR) or by imposing weak disposability and directional distances (nonparaeff, Benchmarking). Neither approach enforces the materials balance as a constraint. pgt does, and puts the competing axiom systems behind one interface.

Building a technology

The constructor bundles the data with the flow coefficients. The shipped steeldemo data is a synthetic panel of steel plants on two production routes, with emissions in CO2 units.

data(steeldemo)
tech <- pgt_tech(
  x = steeldemo[, c("coal_coke", "other_fuel", "raw_material", "flux")],
  y = steeldemo$production,
  b = steeldemo$emissions,
  v = 0.01467,
  group = steeldemo$route,
  id = steeldemo$plant
)
tech
#> Pollution-generating technology
#>   DMUs: 180   inputs: 4   good outputs: 1   bad outputs: 1
#>   inputs: coal_coke, other_fuel, raw_material, flux
#>   material flow coefficients: u = [1, 1, 1, 1], v = 0.01467
#>   groups: Integrated (105), Minimill (75)
#>   materials balance: satisfied for all DMUs

The inputs already sit in pollutant-potential units, so the default \(u = 1\) applies to each; v = 0.01467 is the carbon retained per unit of crude steel (the synthetic generator’s assumption; see ?steeldemo for the derivation).

Auditing the identity before estimation

The weak-G-disposability model gives every DMU a self-reference (its own observation as a feasible solution to its programme) only when its own materials balance holds. mb_check() reports the closure gap per DMU and flags the accounts where more pollutant leaves than enters.

mb <- mb_check(tech)
head(as.data.frame(mb))
#>     id      group potential retained       b       gap    rel_gap violated
#> 1 P001 Integrated   4856127 29930.40 4431848 394348.99 0.08120648    FALSE
#> 2 P001 Integrated   4933440 32542.19 4810748  90149.86 0.01827322    FALSE
#> 3 P001 Integrated   4920386 30214.73 4547577 342594.99 0.06962766    FALSE
#> 4 P002 Integrated   3846987 29343.62 3452634 365009.38 0.09488189    FALSE
#> 5 P002 Integrated   3470487 25188.10 3206835 238464.15 0.06871201    FALSE
#> 6 P002 Integrated   3362598 25968.38 3162682 173948.43 0.05173036    FALSE
attr(mb, "n_violations")
#> [1] 0

An infeasible weak-G-disposability program always belongs to a DMU with a strictly negative gap, so this audit locates any estimation problem in advance. Such DMUs return NA scores and pgt() warns; all other DMUs score in \((0, 1]\).

Estimating environmental efficiency

The weak-G-disposability model minimises the bad output subject to output and input feasibility, the peer emission envelope, and the DMU’s own materials-balance cap (the programmes are stated in full in vignette("models", "pgt")). The score is \(b^*_l / b_l\).

fit <- pgt(tech, model = "wgd")
summary(fit)
#> pgt fit: model = wgd, returns = vrs, peers = all
#>   pollutant: b
#>   DMUs: 180   failed LPs: 0
#> 
#> Efficiency (b*/b):
#>     0%    25%    50%    75%   100% 
#> 0.1269 0.1639 0.4850 0.7122 1.0000 
#> 
#> By group:
#>       group   n n_na median   mean median_dual_output
#>  Integrated 105    0 0.1784 0.3110             0.4733
#>    Minimill  75    0 0.6273 0.7024             0.3018

The envelope model frees the inputs and reduces the program to the convex lower envelope of the \((y, b)\) scatter. Its self-reference is always feasible.

fit_env <- pgt(tech, model = "envelope")
summary(fit_env)
#> pgt fit: model = envelope, returns = vrs, peers = all
#>   pollutant: b
#>   DMUs: 180   failed LPs: 0
#> 
#> Efficiency (b*/b):
#>     0%    25%    50%    75%   100% 
#> 0.1269 0.1639 0.4850 0.7122 1.0000 
#> 
#> By group:
#>       group   n n_na median   mean median_dual_output
#>  Integrated 105    0 0.1784 0.3086             0.4733
#>    Minimill  75    0 0.6273 0.6871             0.3018

Shadow prices and abatement cost

The output-constraint dual \(\eta_l = \partial b^*/\partial y\) (returned as dual_output) measures the marginal emission content of the good output along the frontier. mac_curve() turns it into a marginal abatement cost curve: reducing the bad output by one unit costs \(1 / \eta_l\) units of the good output, or \(p / \eta_l\) at output price \(p\). Read the curve with care about which margin is which: each DMU’s plotted abatement is its distance to the frontier, \(b_l - b^*_l\), which the model itself prices at zero output loss, while the mac value is the marginal cost of abating beyond the DMU’s frontier point. The curve therefore ranks DMUs by the shadow price at their projection; the area under it is not a total-cost estimate.

head(shadow_prices(fit))
#>     id      group       b   b_star dual_output mb_headroom
#> 1 P001 Integrated 4431848 590172.6   0.3017576     4236024
#> 2 P001 Integrated 4810748 643896.2   0.3017576     4257002
#> 3 P001 Integrated 4547577 596021.2   0.3017576     4294150
#> 4 P002 Integrated 3452634 578102.7   0.3017576     3239540
#> 5 P002 Integrated 3206835 492624.8   0.3017576     2952674
#> 6 P002 Integrated 3162682 508675.1   0.3017576     2827955
mac <- mac_curve(fit, price = 550)
plot(mac)

Step curve of marginal abatement cost against cumulative abatement potential

Decomposing across production routes

With a technology group, pgt_decompose() splits environmental efficiency into a within-group and a technology-gap component. The envelope decomposition is an exact identity, Total \(=\) WR \(\times\) TGR, where WR is within-route reallocation and TGR the technology-gap ratio of the metafrontier tradition (Battese, Rao and O’Donnell 2004; O’Donnell, Rao and Battese 2008).

dec <- pgt_decompose(tech, type = "envelope")
summary(dec)
#> pgt decomposition: type = envelope, returns = vrs, DMUs = 180
#>   components: WR x TGR (product = total)
#> 
#> Group medians (components, total, attribution shares):
#>       group   n n_na    WR  TGR total share_WR share_TGR
#>  Integrated 105    0 0.804 0.21 0.178    0.285     0.715
#>    Minimill  75    0 0.627 1.00 0.627    1.000     0.000
#> 
#> (share_* are medians of DMU-level shares; they need not sum to 1.
#>  n_na counts DMUs with any NA component.)

A worked example on real data

The frontier package ships a balanced panel of Philippine rice farms with a physical fertiliser input (NPK, kilograms of active ingredient) and a physical output (PROD, tonnes of paddy). Treating the nitrogen carried by the fertiliser as the pollutant potential gives a nitrogen materials balance. The nitrogen fraction below is illustrative (urea-equivalent, about 0.46 of the active ingredient).

data("riceProdPhil", package = "frontier")
d8 <- riceProdPhil[riceProdPhil$YEARDUM == 8, ]
uN <- 0.46
rice <- pgt_tech(
  x = as.matrix(d8[, c("AREA", "LABOR", "NPK", "OTHER")]),
  y = d8$PROD,
  b = uN * d8$NPK,
  u = c(AREA = 0, LABOR = 0, NPK = uN, OTHER = 0),
  v = 0,
  id = as.character(d8$FMERCODE)
)
rice_fit <- pgt(rice, model = "wgd")
summary(rice_fit)
#> pgt fit: model = wgd, returns = vrs, peers = all
#>   pollutant: b
#>   DMUs: 43   failed LPs: 0
#> 
#> Efficiency (b*/b):
#>     0%    25%    50%    75%   100% 
#> 0.1341 0.4102 0.5939 1.0000 1.0000

Here environmental efficiency is the scope to cut fertiliser nitrogen while holding rice output at its observed level. One caveat on this construction: because the bad output is defined as the full nitrogen potential of a single input (\(b = u_N \cdot\) NPK), every farm’s account closes with equality by construction, the materials-balance cap can never bind, and mb_check() passes trivially; the example illustrates the workflow, not the audit.

The other vignettes state the estimating programmes in full (vignette("models", "pgt")), reproduce the published-result replications (vignette("replication", "pgt")), compare the axiom systems (vignette("comparing-axioms", "pgt")), demonstrate multi-pollutant technologies and heterogeneous coefficients (vignette("multiple-pollutants", "pgt")) and cover productivity measurement with inference (vignette("productivity", "pgt")).

References