Every LP kernel now checks the point lp_solve returns against the
rows and bounds of its programme. lp_solve accepts a point once each
row holds within its tolerances on its own internally scaled model,
and on the unit-magnitude copy of the data a small unit's rows are of
order 1e-7, where those tolerances can exceed several per cent of
the row. On uscoal with the sulfur coefficients scaled by 0.9,
model = "wgd_input_fixed" returned status 0 at Painesville (plant
2936, 38 MWh) with b_star 6 per cent above the cap (mb_headroom
-0.39 t), although no peer mix emits less than the plant's own
6.543 t and the programme is infeasible. A row now fails when its
residual exceeds 1e-5 times the largest of its right-hand side, its
largest term and the unit's own pollutant level, plus 1e-7 times the
total peer weight times the row's largest coefficient on a peer
weight in the solution (lp_solve leaves weights of about 1e-8 on
peers it does not use); a failing point is reported as infeasible
(status 2) with NA scores, and a unit whose own observation is
feasible is retried under the alternative lp_solve settings as
before. On the shipped data no other status or score changes: across
every model, both peer sets and returns to scale, the decompositions
and 50 perturbations of the uscoal sulfur coefficients, the only
change is Painesville's status in the 0.9 scenario and in the 13
draws where it was reported above its cap.
model = "wgd_input_fixed" now decides whether a unit meets its
materials-balance caps (and so whether an infeasibility verdict is
retried) on the data with the relative 1e-5. The previous test
allowed an absolute 1e-8 on the unit-magnitude copy, 0.78 t of SO2 on
uscoal, so it took Painesville's 0.39 t excess for a met cap,
rejected lp_solve's correct verdict of infeasibility and accepted the
point above the cap from the first alternative scaling.
Documentation: citations corrected after a check against the cited papers. No estimator changes.
"wgd" scores it by the minimal-emission programme
of Rodseth (2025, Eq. 6).model = "fdmo" is Rodseth's Eq. 13 with abatement fixed at the
DMU's own level; Eq. 14 is the equivalent convex-regression form
(Proposition 4). Version 0.7.0 cited Eq. 14 for the programme.pigfarms replication (vignette and tests) uses the good-output
coefficient v = 7/6 implied by farms A and B, under which "wgd"
reproduces every Table 2 target of Rodseth (2025). The v = 0 mapping
left farms A and B tied as peers, so the targets of farms B and C
were solver-dependent; that was a property of the mapping, not of
the source.mean_eff of "byprod" is described as Murty, Russell and
Levkoff's by-production efficiency index (one good and one bad
output), not as a variant of the graph measure.boot_pgt(): Kneip, Simar and Wilson (2008) prove the m-out-of-n
bootstrap with replacement consistent under VRS; the default rates
are credited to Kneip, Park and Simar (1998) for VRS and to Park,
Jeong and Simar (2010) for CRS.pgt_ml(): the global frontier removes the cross-period
infeasibility (Oh 2010), not the inconsistency of Aparicio, Pastor
and Zofio (2013); BPC is the best-practice gap change; the default
"wd" technology reproduces Oh's constant-returns technology with
returns = "crs".pgt_decompose(): TGR is the technology gap ratio of Battese, Rao
and O'Donnell (2004), the metatechnology ratio of O'Donnell, Rao and
Battese (2008); WR is the efficiency against the group frontier."wgd" fits the
technology of Hampf and Rodseth (2015) with the minimal-emission
programme of Rodseth (2025, Eq. 6), whereas Hampf and Rodseth score
by a ratio of output to emission. For "fdmo", Eq. 13 is the
indirect form of the directional programme of Eq. 12 under common
coefficients, and the sign restrictions come from Eq. 14.
"byprod" cites Murty, Russell and Levkoff (2012, Eqs. 4.5, 4.7,
5.7 and 7.1): their technology has constant returns, and their
application omits the abatement output, as "byprod" does.
"wgd_input_fixed" is documented as differing from weak
G-disposability at fixed inputs. The uscoal help notes that about
95 per cent of fuel sulfur converts to SO2 (Rodseth 2017). The
boot_pgt() and pgt_ml() help state the conditions of the cited
results.Numerical robustness: every solve of "wgd", "envelope",
"wgd_input_fixed" and the stages of pgt_decompose(type = "rodseth") is now checked against the bound that self-reference
gives. When the evaluated unit is among its peers (and, for
"wgd_input_fixed", meets its own materials-balance caps), its
observation is feasible and attains b, so a solve returning
b* > b, or infeasibility, is a suboptimal vertex or a failed solve
of lp_solve. Such a solve, and any numerical failure (status 5), is
rebuilt and re-solved under two alternative lp_solve scalings and
then the primal simplex, a retry the decomposition stages alone had
before and only for failed solves. A unit whose every attempt still
returns b* > b gets status 1 (lp_solve's SUBOPTIMAL code) and NA
scores, and pgt() counts such units in its warning. A solve that
passes the check at the first attempt is unchanged. On a
1,142-plant-year steel panel in tonnes, 0.7.0 reported up to three
infeasible "wgd_input_fixed" programmes per specification as
numerical failures (status 5); every infeasible programme now
reports status 2, and every status matches an independent GLPK
solve.
New test-scaling.R on steeldemo with every quantity multiplied
by 1e6 (up to 1e13). It checks the decomposition's first stage
against its reduced form, solved separately on a copy in million
tonnes; "wgd", "wgd_input_fixed", "envelope" and the
decomposition against a common rescaling of the quantities;
"wgd_input_fixed" feasibility with and without
materials-balance violators; and the retry and its status-1
report. Without the rescaling of 0.7.0, 96 of the 180 first-stage
optima are suboptimal at this magnitude, two with production
technical efficiency above 1, and three self-feasible
"wgd_input_fixed" programmes are reported infeasible. The tests
mock an internal kernel, so testthat (>= 3.1.7) is suggested.
model = "fdmo" adds the row theta_b <= b (projected emission
non-negative). Rodseth's Eq. 14 is derived for coefficients common
to all units; with producer-specific coefficients the balance row
applies the evaluated unit's coefficients to a peer mix of units with
other coefficients, and a unit whose inputs carry little pollutant
could borrow a pollutant-rich peer's output and receive a negative
projected emission (two closed accounts suffice). Under common
coefficients and closed accounts the row never binds, so the
pigfarms Table 3 replication and every closed-account result with
common coefficients are unchanged.
model = "wgd": z_star now adds the retained content of the output
overshoot, z_star = sum(lambda * z) + v_i'(sum(lambda * y) - y_i),
so that z_star - a_star = b_star whenever the peers' accounts
close. The earlier z_star was the peer mix of caps alone and did
not match b_star when the optimal mix produced more than the unit
with v > 0. b_star and the scores are unchanged, and so is the
pigfarms Table 2 replication (v = 0).
compare_models() ranks scores within 1e-8 of each other as ties
in the Spearman matrix and the bottom-quartile overlap. Efficient
units score 1 up to solver noise (below 1e-9 on the shipped data),
and ranking that noise made the correlations depend on the solver
build: mb_cost-wd on uscoal was 0.221 with one 0.7.0 build and
0.247 with another. The returned scores are unrounded.
The Monte Carlo study inst/simulations/coverage.R now records
interval availability (boot_pgt() returns NA bounds when fewer
than B/2 subsample LPs are feasible, which under VRS hits the units
with the largest output), reports coverage among units with an
interval and the share of all units whose interval exists and
covers, adds a curved frontier to the affine one and an m and
kappa sensitivity at L = 100, keeps per-replication results in
coverage-reps.csv, and reports Monte Carlo standard errors across
replications.
Documentation: the unit dependence of the "fdmo" gross score ends
at the reported sum. The balance row fixes theta_b given
theta_y, so the programme maximises theta_y and a coherent change
of units leaves the optimum unchanged; the earlier help text said
the optimiser could change. The replication vignette no longer calls
the farm E row of Rodseth's Table 3 inconsistent with a common v:
the paper prints one decimal, and common coefficients within that
rounding reproduce every printed row, as a new test asserts.
New dataset with an observed abatement output, referee-driven
documentation and test additions, and a regenerated
steeldemo panel with abatement. No estimator changed, but every
result computed on steeldemo differs from 0.6.1.
steeldemo now carries an observed abatement output and its
pollution-control input: captured (CO2 leaving through
post-combustion CCS at eight integrated plants adopting in 2022 or
2023, capture for utilisation at six plants, or slag mineralisation
at ten), capture_energy (the CO2 potential of the energy the
capture consumes, the natural x_abate marker) and abatement_tech.
The materials-balance identity now closes exactly in every row,
u'x - v y = emissions + captured, so the directional model and the
abatement stages of the five-component decomposition can be run on
the panel. The generator's seed and plant structure are unchanged,
but the earlier process-recovery slack is gone and the emissions
column is re-drawn, so scores computed on steeldemo with 0.6.1 do
not reproduce on this version (the shipped uscoal and pigfarms
data are untouched).
New dataset usfgd: 154 US coal-fired power plants with flue-gas
desulphurisation in 2023, built from EIA-923 (fuel, sulfur,
generation, Schedule 8C control units, sorbent, FGD electricity) and
EPA CAMD monitored SO2, with the SO2 removed by the scrubbers as a
derived abatement output (potential minus measured SO2) and the
reported removal efficiency as a cross-check. The account closes by
construction, so the directional model and the abatement stages of
the five-component decomposition run on real data.
Numerical robustness: pgt(), boot_pgt(), pgt_decompose() and
pgt_ml() now solve every programme on a copy of the technology
scaled to unit magnitude and scale level-valued results back. The
programmes are homogeneous of degree one in the quantities, so
weights, ratio scores and duals are unchanged in exact arithmetic;
at tonne magnitudes lp_solve's default scaling could fail, or stop
at a suboptimal vertex with status 0 (the decomposition's first
stage returned production technical efficiency above 1 for three
steeldemo plants built without the abatement output), and the
rescaled solves remove both. Results can differ from 0.6.1 at the
solver's tolerance.
New test-invariance.R: DMU-order invariance and duplicate-DMU
neutrality for the wgd, wgd_input_fixed and envelope models,
and a two-unit example pinning the retained-content dichotomy (a
unit's wgd score either coincides with its envelope score for
every v or moves for every positive v, switching peers once v
exceeds the emission boundary's descent rate).
?mb_check documents the zero-potential case: rel_gap is computed
against a machine-epsilon floor rather than returned as -Inf, and
such an account is flagged whenever b exceeds tol times machine
epsilon.
?pgt states explicitly that the by-production mean_eff column is
a package summary, not a measure defined in Murty, Russell and
Levkoff (2012).
Documentation: the directional model "fdmo" is now cited as
Rodseth (2025) Eq. 14, the factorially determined multi-output form
that his Proposition 4 derives from Eq. 13 by fixing abatement at the
DMU's own level. Earlier versions cited Eq. 13 for the same programme;
no code changed.
Documentation: DESCRIPTION, the README, ?pgt and the vignettes no
longer describe the materials-balance account as "enforced". The
account enters each programme differently (a cap in
"wgd_input_fixed", an equality in "fdmo", a price system in
"mb_cost", an objective valuation only in "wgd" and
"envelope"), and mb_check() is where it is audited. No code
changed.
Documentation scoping and validation hardening after a full paper-code audit; no estimator changes.
mb_check()'s count of
accounts with a strictly negative closure gap is renamed from
n_exact to n_negative — the old name read as a count of exactly
closing accounts, the opposite of what it counts. The
tolerance-flagged count n_violations is unchanged.?uscoal documents the SO2 measurement provenance (eGRID values
are stack measurements only for CAMD-reporting units, and
emission-factor estimates built from EIA heat input elsewhere, with
a CHP plant's SO2 allocated to electricity) and the total-fuel
basis of coal, which at combined-heat-and-power plants tilts
accounts toward looser closure.pgt() documents the directional model's open-account warning
threshold (relative closure gap above 1e-6), mb_cost's status 2
without a solver call for units with no pollutant-bearing inputs,
and that "wd" reports no output duals; the documented
boot_pgt() default m gains its max(2, .) floor.wgd convention for the implied uncontrolled
emission on deliberately open accounts: z_star is the peer mix of
accounting caps sum lambda (u'x - v'y), which equals
sum lambda (b + a) only when the peers' accounts close.Weak-G-disposability in its Eq. 6 reduced form, several intended outputs, and the five-component decomposition of Rodseth (2025), from a close reading of the published paper.
pgt(model = "wgd") now implements Rodseth (2025) Eq. 6/7 as
printed: only the intended outputs are held at the evaluated unit's
levels; inputs are decision variables. The reduced form minimises
sum lambda (b_l + v_i'y_l) - v_i'y_i over output-constrained peer
mixes. The programme is always self-feasible, scores lie in (0, 1]
on any data, and materials-balance violations no longer affect it
(they remain a data audit in mb_check()). When an abatement output
is observed the fit reports the projection's implied uncontrolled
emission z_star and abatement a_star, the quantities of the
paper's Table 2, which the test suite now replicates. Output duals
are total derivatives and can be negative under large retained
content. Scores on any data with v > 0 differ from 0.5.x.model = "wgd_input_fixed" (briefly named
wgd_anchored in unreleased development commits), documented as the
package's own current-input benchmark rather than Eq. 6. It keeps
the cap, mb_headroom, infeasibility confined to violators, and the
per-pollutant cap rows.pgt_tech() accepts y as an L x M matrix with named columns and
v in per-output shapes up to an L x M x P array; every historical
single-output shape still works. All estimators except fdmo and
pgt_ml handle M > 1 (those two stop with a clear error, since
their directions are defined for a single output). Fits carry one
column per output and per-output duals; mac_curve() requires a
single-output fit.pgt_decompose(type = "rodseth") is now the five-component
decomposition of Rodseth (2025, Eq. 11), built on stage programmes
of the extended representation (Eq. 9): TE production, Quality
(producer-specific coefficients), AE production, TE abatement
(observed abatement output) and AE abatement (dedicated
pollution-control inputs). Components collapse to exactly 1 when the
data cannot separate them, and their product telescopes to the wgd
efficiency. Technology groups are no longer required for this type.
The stage programmes use measurement-consistent uncontrolled
emissions z = b + a, which coincide with the paper's accounting
when accounts close.x_abate is now functional: it partitions the inputs for the two
abatement stages, and the recorded-but-unused warning is gone.lp_solve scalings, each attempt
time-limited: with the default scaling alone, self-feasible stages
on large-magnitude data (emissions of order 1e7) were reported as
numerically failed for a handful of units.pgt_ml() now defaults to technology = "wd", the
weak-disposability reference under which Oh (2010) defines the
index, and the input-keeping exploratory option is named
"input_envelope" (previously "wgd"): it keeps the input rows of
the pre-0.6.0 programme, so after the wgd redefinition its old label
no longer matched the estimating model of the same name. No
technology label is now shared with a pgt() model whose programme
differs.inst/simulations/coverage.R writes its results next to the script
when the package-source path is absent, so a sourced installed copy
no longer fails at the final write.boot_pgt() returns NA bounds and standard error for a DMU whose
feasible replicates fall below half of B, instead of an interval
resting on the remainder; n_ok still reports the count.print.pgt_tech() separates the dimension counts with semicolons
and wraps the DMU-specific-coefficient and audit lines within a
76-character console.b_star warning of model = "mb_cost" no longer
attributes the sign to DMU-specific coefficients, which need not be
present; any positive retained-content coefficient permits it.uscoal: 212 US coal-fired power plants (2022) with a
measured sulfur-dioxide materials-balance account, built from
EIA-923 and EPA eGRID2022 (construction in data-raw/uscoal.R).
Producer-specific sulfur-based coefficients, an FGD sorbent
pollution-control input for x_abate, an FGD/no-FGD group, and 8
genuinely violating accounts for the audit. ?uscoal documents that
sorbent is 0 both without SO2 controls and at the 36 FGD plants
that report no sorbent quantity, so fgd, not a zero sorbent, is
the scrubbing indicator.plot.pgt_compare() varies line type as well as colour, so the
series stay readable in greyscale, and labels the y axis
"Principal score (b*/b; EE for mb_cost)" when mb_cost is among
the compared models: its EE is a material-inflow ratio, not an
emission ratio. Passing lty, col, xlab or ylab overrides
the defaults.mean_eff (previously fgl): the column is a plain average,
not the Fare-Grosskopf-Lovell hyperbolic graph measure, and its
old name invited the confusion.status column is documented as the lp_solve
solver code (0 solved, 2 infeasible, 5 numerically failed), and the
directional model's documentation warns that rescaling y or b
changes its implicit direction weighting, not merely the reported
number.efficiency column is called
the principal environmental-efficiency score throughout the
documentation and printed output; the internal helper was renamed
accordingly.Revisions from a full peer review of the package as a software submission.
boot_pgt() no longer forces the evaluated DMU into every subsample:
each unit is scored against the subsample frontier alone, so frontier
units carry genuine resampling variation instead of a degenerate
zero-width interval. Replicates with an infeasible programme surface
as NA; the new per_dmu$n_ok column counts feasible replicates and
a warning fires when they run thin. Interval endpoints truncate to
[0, max(1, estimate)], respecting above-1 violator scores.kappa for model = "byprod" now uses the
polluting-input dimension of the T2 sub-programme; boot_pgt() warns
when a panel technology is subsampled as a cross-section; group-mean
intervals are documented as a descriptive stability band (the
aggregate theory of Kneip, Simar and Wilson 2015 is not implemented).boot_pgt_sensitivity() reports the implied rate exponent
kappa_hat from the Politis-Romano-Wolf log-spread regression, so
the borrowed DEA rate is checkable, plus the median over strictly
positive widths.inst/simulations/coverage.R; its results are reported in the
productivity vignette.pgt_ml() gains technology = "wd": the weak-disposability
directional distance in the Kuosmanen form, the technology under
which Chung, Fare and Grosskopf (1997) and Oh (2010) define the
index, making it the comparator defined by Oh (2010). The free-disposal
semantics of the default "wgd" and "envelope" options are now
stated plainly and the claim that they coincide with the
weak-disposability index when the cap does not bind has been removed
(it was not what the kernel computed).models states every estimating programme in
full, and multiple-pollutants demonstrates multi-pollutant
technologies, DMU-specific coefficients and the staged decomposition.compare_models() documents that the models' principal scores are
different
quantities (EE for mb_cost) and that only rank-based statistics are
strictly comparable; mac_curve() documents which margin the curve
prices; foundational and origin literature is cited (Ayres and
Kneese 1969; Lauwers 2009; Dakpo, Jeanneaux and Latruffe 2016; Chung,
Fare and Grosskopf 1997; Chambers, Chung and Fare 1996; Kuosmanen and
Podinovski 2009; O'Donnell, Rao and Battese 2008; Hampf 2014; Fare,
Grosskopf, Lovell and Yaisawarng 1993), every vignette carries a
reference list, and DOIs are attached throughout.inst/CITATION added; a package overview with a scope statement
replaces the empty pgt-package help page; the README states the
package's scope and labels the experimental and heuristic modules.Fixes from a whole-package review.
boot_pgt() now builds proper subsampling intervals (Politis and
Romano 1994): the subsample distribution is recentred at the point
estimate and rescaled by the relative convergence rate (m/L)^kappa,
so intervals extend downward from the estimate, matching the
direction of the frontier bias. Previous versions built percentile
intervals from the raw subsample distribution, which lies entirely at
or above the point estimate. The reported se is rescaled the same
way. A new kappa argument controls the rate exponent, defaulting to
the DEA rate for the model's effective dimension.boot_pgt() restores the caller's RNG state on exit instead of
permanently reseeding the session.pgt(model = "wgd") now enforces every pollutant's
materials-balance cap on the peer mix, as documented; earlier
versions constrained only the selected pollutant.pgt(model = "fdmo") warns when material accounts do not close
exactly (u'x - v y != b + a): the model imposes the identity as an
equality, so open accounts either make the LP infeasible or shift the
closure gap into the scores. mb_check() reports the equality
residual in a new closure column when abatement is observed.pgt(model = "byprod") or "mb_cost" now sets
every score of that DMU to NA, so status != 0 always means the
scores are NA, matching the warning text and boot_pgt().pgt(model = "mb_cost") warns when the implied minimal emission
b_star is negative (possible under DMU-specific coefficients).pgt() and pgt_decompose() warn that x_abate is recorded but not
yet used by any estimator.pgt_tech() rejects an ambiguous v whose length equals both the
number of DMUs and the number of pollutants, and rejects duplicated
(id, period) pairs, which pgt_ml() silently mismatched before.compare_models() recognises the documented aliases: wgd_rodseth
is accepted, and ddf now yields the clear directional-model
rejection instead of a misleading unknown-model error. Duplicated
entries are ignored with a warning, and the largest-disagreement line
no longer reports a model's self-correlation.mac_curve() and shadow_prices() now error clearly for models
without output duals (fdmo, byprod, mb_cost, wd); earlier
versions returned an empty curve reporting zero exclusions.print() on a subset of an mb_check() audit reports counts for the
printed rows, not the full table.summary() of a pgt_ml fit no longer errors when transitions
contain NA; failed distance LPs are counted and warned about.envelope) and by summary() as well.pgt(pollutant = <character vector>) gives a clear error instead of
an internal condition-length failure.boot_pgt()).Benchmarking dropped from Suggests (never used).Panel productivity measurement and inference.
pgt_ml(): the global Malmquist-Luenberger productivity index (Oh
2010) under the pollution-generating technology, decomposed into
efficiency change and best-practice change (GML = EC x BPC). The
global reference frontier removes the cross-period infeasibility of the
adjacent-period index. Marked experimental: the index under a
materials-balance technology is not yet settled in the literature.boot_pgt() and boot_pgt_sensitivity(): subsampling (m out of n)
confidence intervals for per-DMU scores and group means, with a
subsample-size sensitivity check. The intervals are heuristic; the
documentation states why.pgt_tech() gains a period argument for panel data.pgt_ml and pgt_boot.Competing axiom systems and a comparison harness.
pgt(model = "byprod"): the by-production intersection technology of
Murty, Russell and Levkoff (2012). Reports output efficiency, emission
efficiency and their graph average. Unit tests reproduce the analytic
scores of the paper's Example 1.pgt(model = "mb_cost"): the materials-balance cost model of Coelli,
Lauwers and Van Huylenbroeck (2007), with the EE = TE x EAE
decomposition into technical and environmental-allocative efficiency.pgt(model = "wd"): the weak-disposability model (Kuosmanen 2005
correct VRS formulation) as a reference axiom system.compare_models(): fits the efficiency-scored models on one technology
and reports Spearman rank agreement, worst-quartile overlap and the
largest ranking disagreement.pgt_tech() gains a polluting argument marking the emission-causing
inputs of the by-production model.The directional representation and richer material accounting.
pgt(model = "fdmo"): the factorially determined multi-output
(directional) representation of Rodseth (2025), Eq. 14, that is
Eq. 13 with abatement fixed at each DMU's own level (alias "ddf"). Jointly expands the good
output and contracts the bad output; returns good-output, bad-output
and gross efficiency. Unit tests reproduce the paper's Table 3.u may be an L x N matrix
and v a length-L vector, so material or output quality can vary
across units (Eder 2022, Rodseth 2025).b may be an L x P matrix, with u a
named list and v a vector or matrix; each pollutant carries its own
materials-balance identity, and a pollutant argument selects the one
to estimate.pgt_tech() gains an a (abatement output) argument, required by the
directional model.Initial release: the Rodseth (2025) weak-G-disposability slice.
pgt_tech(): pollution-generating technology constructor with
materials-balance metadata (material flow coefficients u, v).mb_check(): pre-estimation audit of the materials-balance identity
u'x - v y >= b.pgt(): per-DMU minimum-bad-output estimation; models "wgd"
(Rodseth 2025, Eq. 6 reduced form, alias "wgd_rodseth") and
"envelope" (Eq. 6 with inputs free: the convex lower (y, b)
envelope); VRS/CRS; pooled or group reference sets.pgt_decompose(): metafrontier decompositions of environmental
efficiency; "envelope" (Total = WR x TGR, exact identity, no
feasibility screen) and "rodseth" (TE x Technology x AE from staged
relaxations of the full WGD program).shadow_prices() and mac_curve(): constraint duals and marginal
abatement cost curves.print, summary, plot, as.data.frame methods; synthetic
steeldemo panel; testthat suite with hand-computed LP checks and a
cross-check against an independent reference implementation of the
reduced-form program.pigfarms data: the numerical example of Rodseth (2025), Table 1.
Unit tests replicate the paper's Table 2 minimal controlled
emissions (16, 16, 16, 20, 16) and recover its abatement column from
the mb_check() closure gap.