Multiple pollutants and heterogeneous coefficients

A technology can carry several pollutants, each with its own materials-balance account, and the material flow coefficients can vary across units (heterogeneous input quality; Eder 2022, Rødseth 2025). This vignette demonstrates multi-pollutant accounts, DMU-specific coefficients, several intended outputs and the five-component decomposition, on a small synthetic example.

A two-pollutant technology

Forty plants burn coal and gas. Each fuel carries CO2 potential and SO2 potential, with different coefficients; a share of the CO2 potential is retained in the product (\(v_{CO2} = 0.1\)), none of the SO2. b is a matrix with one named column per pollutant, u a named list with one element per pollutant, and v a named vector.

set.seed(42)
L <- 40
coal <- runif(L, 20, 60)
gas <- runif(L, 10, 40)
y <- runif(L, 15, 45)
u_co2 <- c(coal = 2.5, gas = 1.4)
u_so2 <- c(coal = 0.020, gas = 0.001)
cap_co2 <- 2.5 * coal + 1.4 * gas - 0.1 * y
cap_so2 <- 0.020 * coal + 0.001 * gas
b_co2 <- cap_co2 * runif(L, 0.55, 0.95)
b_so2 <- cap_so2 * runif(L, 0.55, 0.95)
# three plants under-report scrubbing and breach their SO2 account
b_so2[1:3] <- cap_so2[1:3] * runif(3, 1.05, 1.25)

tech2 <- pgt_tech(
  x = cbind(coal = coal, gas = gas),
  y = y,
  b = cbind(co2 = b_co2, so2 = b_so2),
  u = list(co2 = u_co2, so2 = u_so2),
  v = c(co2 = 0.1, so2 = 0),
  group = factor(rep(c("A", "B"), each = L / 2))
)
tech2
#> Pollution-generating technology
#>   DMUs: 40;  inputs: 2;  good outputs: 1;  bad outputs: 2
#>   inputs: coal, gas
#>   material flow coefficients [co2]: u = [2.5, 1.4], v = [0.1]
#>   material flow coefficients [so2]: u = [0.020, 0.001], v = [0]
#>   groups: A (20), B (20)
#>   materials balance: 3 of 80 accounts violate the identity
#>     (see mb_check())

mb_check() audits one account per plant and pollutant; the three seeded SO2 violations surface immediately.

mb <- mb_check(tech2)
table(subset(as.data.frame(mb), violated)$pollutant)
#> 
#> so2 
#>   3

Every pollutant’s cap constrains the input-fixed benchmark

The faithful "wgd" programme of Eq. 6 uses only the selected pollutant’s data. Cross-account discipline lives in the input-fixed benchmark: pgt(model = "wgd_input_fixed") selects one pollutant for the objective through pollutant, but the caps of every pollutant constrain the peer mix (see vignette("models", "pgt") for the collapsed cap rows). A peer that looks attractive on CO2 can be unusable because leaning on it would breach the evaluated plant’s SO2 account.

fit_co2 <- pgt(tech2, model = "wgd_input_fixed", pollutant = "co2")
summary(fit_co2)
#> pgt fit: model = wgd_input_fixed, returns = vrs, peers = all
#>   pollutant: co2
#>   DMUs: 40   failed LPs: 0
#> 
#> Efficiency (b*/b):
#>     0%    25%    50%    75%   100% 
#> 0.2716 0.3897 0.5354 0.8206 1.0000 
#> 
#> By group:
#>  group  n n_na median   mean median_dual_output
#>      A 20    0 0.4644 0.5236              2.645
#>      B 20    0 0.6411 0.6751              2.716

Dropping the SO2 account (a single-pollutant technology on the same data) weakly lowers every b_star, since the projection loses one set of constraints:

tech1 <- pgt_tech(
  x = cbind(coal = coal, gas = gas), y = y, b = b_co2, u = u_co2,
  v = 0.1, group = factor(rep(c("A", "B"), each = L / 2))
)
fit1 <- pgt(tech1, model = "wgd_input_fixed")
summary(fit1$results$b_star - fit_co2$results$b_star)
#>       Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
#> -1.421e-14  2.842e-14  4.263e-14  8.148e-13  6.573e-14  1.027e-11

The three SO2 violators show the documented boundary behaviour: a plant violating another pollutant’s identity loses its self-reference, so its CO2 score can exceed 1 (its MB-consistent projection emits more CO2 than the plant reports) or its programme can be infeasible. Neither happens for plants whose accounts all hold.

head(cbind(fit_co2$results[c("id", "efficiency", "status")],
           so2_violated = mb$violated[mb$pollutant == "so2"]), 5)
#>   id efficiency status so2_violated
#> 1  1  0.4613612      0         TRUE
#> 2  2  0.3388767      0         TRUE
#> 3  3  0.8210681      0         TRUE
#> 4  4  0.4116247      0        FALSE
#> 5  5  0.6013043      0        FALSE

DMU-specific coefficients

Coefficients can vary by plant, for example with coal quality (Eder 2022): pass an \(L \times N\) matrix per pollutant instead of a vector. Here the first twenty plants burn higher-carbon coal.

U_co2 <- matrix(u_co2, L, 2, byrow = TRUE,
                dimnames = list(NULL, c("coal", "gas")))
U_co2[1:20, "coal"] <- 2.8
tech_het <- pgt_tech(
  x = cbind(coal = coal, gas = gas), y = y,
  b = pmin(b_co2, (2.8 * coal + 1.4 * gas - 0.1 * y) * 0.95),
  u = U_co2, v = 0.1
)
fit_het <- pgt(tech_het, model = "wgd_input_fixed")
median(fit_het$results$efficiency, na.rm = TRUE)
#> [1] 0.5354456

Each plant’s cap now uses its own coefficients, so two plants with the same fuel bill face different materials-balance ceilings. The input coefficients do not enter the reduced form of the faithful "wgd" programme (only the output coefficients \(v_i\) do); producer-specific \(u\) instead surfaces in the quality component of the decomposition below.

The five-component decomposition

pgt_decompose(type = "rodseth") splits the weak-G-disposability score into the five components of Rødseth (2025, Eq. 11): production technical efficiency, coefficient quality, production allocative efficiency, abatement technical efficiency and abatement allocative efficiency (the stage programmes are stated in vignette("models", "pgt")). A component collapses to exactly 1 when the data cannot separate it. Here no abatement output is observed and no pollution-control input is marked, so the two abatement components are 1, while the heterogeneous coefficients of tech_het activate the quality component.

dec <- pgt_decompose(tech_het, type = "rodseth")
summary(dec)
#> pgt decomposition: type = rodseth, returns = vrs, DMUs = 40
#>   components: te_production x quality x ae_production x te_abatement x ae_abatement (product = total)
#> 
#> Group medians (components, total, attribution shares):
#>  group  n n_na te_production quality ae_production te_abatement ae_abatement
#>    all 40    0         0.991       1          0.57            1            1
#>  total share_te_production share_quality share_ae_production share_te_abatement
#>  0.524               0.127         0.068               0.719                  0
#>  share_ae_abatement
#>                   0
#> 
#> (share_* are medians of DMU-level shares; they need not sum to 1.
#>  n_na counts DMUs with any NA component.)

The product of the five components equals the "wgd" efficiency exactly, and technology groups are not required for this type.

Several intended outputs

Good outputs can be a matrix with one named column per output, and the retained content becomes \(v_i'y_i\) with one coefficient per output. Here the plants sell power and process heat, retaining different carbon shares per unit, and the programme holds both outputs at the evaluated plant’s levels.

y2 <- cbind(power = y, heat = runif(L, 5, 20))
tech_m <- pgt_tech(
  x = cbind(coal = coal, gas = gas),
  y = y2,
  b = b_co2,
  u = u_co2,
  v = c(0.10, 0.05)
)
fit_m <- pgt(tech_m, model = "wgd")
head(fit_m$results, 3)
#>   id    power      heat         b   b_star efficiency dual_power dual_heat
#> 1  1 32.44812 11.641944 116.49452 53.77858  0.4616405  3.5644553  1.291545
#> 2  2 19.73716  7.368202 124.30608 42.33774  0.3405926  0.5122657 -0.050000
#> 3  3 25.77085 11.634870  71.39915 45.21526  0.6332745  0.5122657 -0.050000
#>   status
#> 1      0
#> 2      0
#> 3      0