A sigma_u collapsed to zero was reported as the correct answer even when
the residual skew said inefficiency was present. sfm() sets
$sigma_u_at_bound whenever the fitted one-sided scale reaches the zero
boundary, and the two cases that reach it read in OPPOSITE directions. Under
WRONG skew (third central moment >= 0) the boundary genuinely is the
maximum likelihood estimate -- the Type I failure of Olson, Schmidt and
Waldman (1980) -- and should be read as no evidence of inefficiency. Under
CORRECT skew the data say inefficiency IS present, so the boundary should not
be the answer and the fit is suspect.
Only the first case was handled. The fit-time warning was guarded on
wrong_skew && at_bound, so the second case was flagged and silent.
Worse, print() and summary() printed the sentence "Under wrong skewness
this is the correct MLE, not a failure (Waldman 1982)" in both cases,
having only dropped the clause that names the skew -- so the suspect fit was
shown the exculpation belonging to the sound one, in the output that survives
saveRDS() and a fresh session.
Measured on sfm(model_name = "NW"), y_pcs_w, N = 150: seeds 1, 5 and 7
all return optim() convergence code 52
(ABNORMAL_TERMINATION_IN_LNSRCH), no standard errors,
$sigma_u_at_bound = TRUE, $wrong_skew = FALSE and residual third moments
of -1585.5, -469.3 and -4007.9. Two of the three are not even
stationary: moving sigma_v alone, with every other parameter held at its
fitted value, gains 0.4771 and 0.0452 in log-likelihood. On seed 1 a
point 6.5743 better exists at the same sigma_u. So a user saw a boundary
flag, no standard errors and a printed note calling it correct, on a fit that
had not converged and was not a maximum.
Both reporting paths now branch on the skew. The warning for the suspect case names the residual moment, reports a non-zero convergence code when there is one, says the fit may not be a maximum, and states that the wrong-skew reading does not apply. The sound case keeps the message it had.
"tHN" is excluded from the new reading and keeps its own: it already warns
that heavy-tailed noise absorbing the whole one-sided component is a known
property of the model rather than a numerical failure, and its row names put
it through this same generic block, so the two would otherwise contradict
each other on one fit.
sfm(model_name = "NTN") reported a log-likelihood above the true one, and
above OLS, as lambda approached its lower bound. The truncated-normal
density was formed as log Phi(aa) - log Phi(bb) with
aa = (mu/lambda - eps*lambda)/sigma and
bb = (mu/sigma)*sqrt(1 + lambda^-2). Both arguments carry the same
mu/lambda, so with mu < 0 each term reached about -1e13 while their
difference stayed O(1): the answer was the rounding error of two enormous
numbers. On one sample the reported value was -370.99 where the true value
is -372.73, which also put it 1.21 ABOVE the OLS log-likelihood -- and
near lambda = 0 the model collapses to a normal regression, so OLS is an
upper bound there and cannot be beaten. The error was granular in exactly the
way the diagnosis predicts: on one fit it came out as a multiple of the
floating-point spacing of aa^2/2, 512, once per observation.
aa^2 - bb^2 is (-2*mu*eps + eps^2*lambda^2 - mu^2)/sigma^2 analytically,
with the mu^2/lambda^2 cancelling exactly, so the difference is now taken
through .log_phi_tilt() -- as "NE", "NGE" and "TSL" already do -- and the
quadratic reduces in closed form to -eps^2*(1 + lambda^2)/(2*sigma^2).
That path is taken only where the cancellation is real, namely both
arguments below -1e3 (the switch .log_phi_tilt() itself uses, below which
the direct form is accurate to 1e-11). Everywhere else -- mu > 0, or
lambda away from its bound -- the pre-existing expression is kept
bit for bit, including its association order. That matters: on the flat
sigma_v ridge where lambda -> Inf is a supremum, lambda is not
identified, and a last-bit change to the likelihood moved the reported
lambda by a factor of 20 at an unchanged log-likelihood. 2413 of 2700 fits
are now bit-identical to before. Verified two independent
ways: against the analytic lambda -> 0 limit, which the fix now approaches
at the expected O(lambda^2) rate down to 2e-16 while the old form turned
around below lambda = 1e-4 and reached 0.63; and, where neither logarithm
underflows, against log(Phi(aa)/Phi(bb)) evaluated in level space.
Reported in #55.
sfm(model_name = "NTN") could return a fit far below the nested "NHN".
At mu = 0 the truncated-normal likelihood IS the half-normal one with the
same (lambda, sigma), and mu = 0 is interior -- start_cs() bounds only
lambda and sigma below -- so a maximum below NHN's is a failure rather
than a limit that is merely approached. It happened on 200 of 2700 fits, by up
to 764. NTN is now checked at the end against a fit made with mu held at 0,
and polished from it when that point is better, exactly as "NG" is checked
against shape 1 (#45) and "NNAK" against m = 0.5 (#30) -- a check at the
end rather than another start, because as a start it moves fits that were
never in trouble. On the same 2700 fits this leaves 13 below NHN, by at most
0.39; all 13 sit at lambda above 1300, which is the flat sigma_v boundary
of the A59 case, where a single polish stops short on the ridge. Reported
in #55.
vcov(type = "sandwich") and vcov(type = "clustered") recommended a fallback
that did not work. When the bread is undefined -- the fit carries no
Hessian, or the Hessian is singular -- both paths advised
type = "bhhh" on the grounds that it needs no Hessian. That is true but not
sufficient: BHHH still has to invert the outer product of the scores, and a
fit that reaches these errors is frequently flat in a parameter, which makes
that singular too. On such a fit the error said type = "bhhh" "is defined
here" while vcov(type = "bhhh") raised "the outer product of gradients is
singular", so a user following the advice met a second error.
The recommendation is now made only when BHHH is actually available, and when
it is not the message says so and names the parameter responsible. Found by
measuring test coverage: .sfa_dead_parameters(), which supplies that name,
had no test coverage at all -- 57 of 57 expressions never executed -- so
the whole diagnostic path for these covariances had never run in the suite.
sfm(model_name = "NTN") reported a collapsed noise component by saying
nothing at all. On some samples NTN has no interior maximum in lambda:
the profile log-likelihood rises monotonically toward a finite limit as
lambda grows without bound, so the fit correctly runs to lambda of order
1e6 and the boundary is the maximum likelihood estimate. But the
likelihood is flat in lambda there, so the Hessian is singular and every
standard error comes back NA. The user was given lambda = 972468, six
NA standard errors, and no message of any kind.
sfm() now warns and sets ntn_sigma_v_at_bound on the fit, saying that the
boundary is the estimate rather than a failure, that no finite standard error
for lambda exists at that point, and that the result should be read as "no
noise component identified in these data". This is the same treatment
tHN has had since 1.2.0 for the mirror case, sigma_u collapsing onto
zero; the other component was simply never covered.
The criterion is sigma_v/sigma < 1e-3, the same scale-free ratio
thn_sigma_u_at_bound uses. Measured over 12 samples it separates the two
regimes by five orders of magnitude -- the boundary fits sit at 4.9e-08 to
1.3e-06 and the interior fits at 0.126 to 0.442 -- so the threshold
falls in empty space rather than being tuned to the data.
marginal_effects() on a truncated-normal fit (sfm(model_name = "NTN")
with uhet/muhet) now reports every covariate of the mu design. The
covariate list was taken from the sigma_u design alone, so with
uhet = ~ z1, muhet = ~ z2 the z2 effect was silently missing, and with
muhet alone the call stopped claiming every column was constant. The
chain rule already handled a covariate in one design only; it was never
given the names. ?marginal_effects also now lists the "component"
attribute and the _h column names returned for component = "h".
sfm(model_name = "NG") could stop at the OLS point, below the NE
solution it nests (issue #45). NG at shape 1 is NE, so its maximum can
never be below NE's. But its starting values only hold E[u] at a moment
anchor while sweeping the shape, and on some samples the optimizer stages
ended with sigma_u on its floor, where the likelihood is OLS's. On the
issue's sample NG returned -359.2835 (OLS's value, every firm fully
efficient) against NE's -359.2456. NG now also fits its likelihood with the
shape held at 1 -- NE's likelihood -- from four splits of the residual
variance, and if the finished fit is below that point it is polished from
there and the better kept, the same check gap A25 added against the
multistart point. On that sample NG now reaches -358.568, confirmed by
numerical integration of the normal-gamma density. Over 360 fits (20
samples of n = 200, every data_gen_cs() column) NG ended more than 1e-3
below NE on 4 and now on none; the other 353 fits are bit-for-bit unchanged.
ivsfm(instruments = ~1) failed with an unrelated error before R 4.4.
With no excluded instruments the instrument matrix was built by
reformulate(character(0)), which R 4.3 and earlier reject ("'termlabels'
must be a character vector of length at least one"). The intended message
-- that the order condition needs at least as many excluded instruments as
endogenous variables -- was reached only on R 4.4 and later. The package
supports R 4.0 and up; the message is now the same on all of them.
ivsfm(model_name = "C2SLS") silently ignored uhet. C2SLS is 2SLS
with a homoskedastic moment correction for sigma_u and the intercept, so
it has no delta to estimate. It accepted uhet anyway and returned the
homoskedastic fit unchanged -- the same coefficients and the same jlms as
without it, and no delta_ rows -- so a user asking for the APS (2017)
model got the 2016 one without being told. Supplying uhet with
"C2SLS" is now an error that points to "IVLIML" or "IVCF".
ttsfm(model_name = "TTNLS") reported wrong standard errors for every
coefficient. Its standard errors came from inverting the Hessian of the
sum of squares directly. That Hessian is 2 X'X in the coefficients, so
the covariance is 2 s^2 times its inverse, with s^2 the residual
variance. Without that factor the slopes' standard errors did not depend on
the noise in the data at all: on 3,000 simulated observations they were
about half of OLS's (0.019 where OLS gives 0.039). The Hessian also kept
the unidentified sigma_u and sigma_w directions, along which the
intercept trades off one for one, so the intercept's standard error depended
only on the sample size (40.8 at n = 300, 12.9 at n = 3,000), whatever the
data. Both are fixed. Without determinants TTNLS is OLS, and its
coefficients and standard errors now equal lm()'s. Where sigma_u and
sigma_w carry determinants the standard errors are conditional on the
fitted determinant coefficients.
Method now rejects any value other than "L-BFGS-B" instead of silently
optimizing without bounds. Method is passed to optim() for the final
stage together with the lower and upper bounds that hold the scale parameters
positive -- but optim() applies bounds for "L-BFGS-B" only. With
"BFGS", "CG", "Nelder-Mead" or "SANN", R warned
("bounds can only be used with method L-BFGS-B") and the stage then ran
UNBOUNDED, free to take a standard deviation to or past zero. That was
invisible on an interior optimum -- all five methods agree to eight figures on
a well-behaved fit, because the bounds never bind there -- and it bit only
where a bound was active, which is exactly the boundary-collapse region the
scaffold's guards exist for. sfm(), zsfm(), lcsfm(), psfm(),
ttsfm(), copsfm(), ivsfm() and selsfm() now validate Method before
fitting anything, and opt.optim() enforces it as well, since that is the
function which passes the bounds.
marginal_effects() now covers the truncated normal. It previously
handled the families whose E[u] is proportional to the scale -- the
half-normal and the exponential -- where the effect collapses to
half * delta_k * E[u]. model_name = "NTN" fitted with uhet/muhet is
not one of those: the pre-truncation mean enters E[u] beside the scale, so
the effect is the chain rule through both designs. With a = mu/sigma_u and
lambda(a) = phi(a)/Phi(a), dE[u]/dmu = 1 - lambda(lambda + a) and
dE[u]/dsigma_u = lambda + a lambda(lambda + a); mu carries an identity
link, so delta_mu is not a semi-elasticity the way the scale coefficients
are. A covariate in one design but not the other contributes only through
that one. Every derivative is checked against numDeriv in the tests rather
than against a second hand derivation.
marginal_effects() also reports the scaling property. Under
sfm(scaling = ~ z) one factor h = exp(z'delta_s) multiplies both the scale
and the pre-truncation mean, so a = mu/sigma_u is constant in z,
E[u] = h E[u*], and the effect is exactly delta_s_k * E[u] whatever the
family -- with delta_s_k itself the semi-elasticity, so dE_u.dz / E_u is
constant by construction. The quantity was already fitted and simply not
reported: the fit stored no scaling block, because n_blocks does not name
one. sfm() now stores s_spec alongside z_spec, v_spec and mu_spec.
vcov() now offers four covariances, not two. type gains
"sandwich" and "clustered", joining "hessian" and "bhhh", and
vcov() gains a cluster argument. Every maximum-likelihood fit in the
package that retains its likelihood through keep_objective = TRUE can
therefore report Hessian, outer-product, heteroskedasticity-robust and
cluster-robust standard errors through one interface, without the caller
having to attach sandwich and know its idiom.
These are computed in the package rather than delegated, because the
covariance has to be formed on the scale the scores live on and only then
mapped to the scale coef() reports. sandwich::sandwich() cannot take
that step -- it composes bread and meat itself and knows nothing about the
mapping. vcov(type = "clustered") clusters over the rows of the score
matrix, which for the panel likelihoods is one row per FIRM rather than per
firm-year, and says so when handed a vector of the wrong length.
OPG (BHHH) standard errors are now available from every entry point whose
likelihood is separable per observation. keep_objective = TRUE gains
zsfm(), lcsfm(), ttsfm() ("TTNE", "TTHN"), copsfm(), ivsfm()
("IVLIML", "IVCF") and selsfm(), joining sfm() and psfm(). Those
fits then support vcov(type = "bhhh"), influence_sfa(), the sandwich
methods, TIC() and vuong(). Because the OPG needs no Hessian at all, it
is defined in exactly the cases where the default path fails -- an
indefinite or non-invertible Hessian.
The estimators that maximise no log-likelihood have no score, and now say so
instead of returning a fit whose score-based tools are silently unavailable:
ttsfm("TTNLS") (nonlinear least squares) and ivsfm("C2SLS") (2SLS with a
corrected intercept) warn that keep_objective = TRUE has no effect, as
psfm() already did for its moment-based and FE estimators. npsfm() is
nonparametric and has no likelihood to retain.
Two caveats are documented rather than papered over. selsfm() is a
two-step estimator, so its OPG conditions on the first-stage probit exactly
as its Hessian already does; and for lcsfm("LCM_CN") with
penalty_c > 0 the contributions are those of the unpenalised likelihood,
since a penalty on the class probabilities is not a per-observation quantity.
psfm()'s closed-form panel likelihoods supply per-firm scores too.
keep_objective = TRUE now reaches "GTRE_FML", "TFE_WMLE", "FD",
"PL80", "PL80_MVTN", "BC92", "K1990" and "K1990modified", joining
the four simulated-ML models. Every psfm() estimator that maximises a
log-likelihood can now produce a score matrix, and the remaining ones
("TFE", "SSFE", "SSRE", "SSCRE", "CSS", "LS", "KSS",
"GTRE_SEQ1", "GTRE_SEQ2") are moment-based, FE/LSDV or sequential and
warn that the argument has no effect. Contributions are per firm, so
sandwich::vcovCL() takes a cluster variable with one entry per firm and
the resulting errors are firm-clustered by construction.
This matters most where the Hessian path fails. "FD"'s Hessian is
routinely indefinite, so three of its five standard errors print as NaN;
the OPG returns all five, and they are large -- which is the honest answer
for parameters the near-zero Hessian eigenvalues said were barely
identified.
NG and NNAK: the likelihood was wrong for residuals far below the frontier, and an NNAK fit could end below the half-normal model it nests (issue #30). Three defects, now fixed.
First, .log_pcf() -- the log parabolic cylinder function both densities
use -- still returned the clipped value of the gsl code it replaced for
arguments below z = -37.42, wrong there by hundreds of log units: at
z = -42.9, log D was 240.9 against a true 461.0. At shape 1/2 NNAK
equals NHN in closed form, yet on one sample NNAK evaluated at NHN's
optimum scored -624.8 against -404.9, all of the gap from one outlying
residual. .log_pcf() is now exact at any argument, which is checked against
the closed form for order -1 and against numerical integration.
Second, for z < 0 both densities formed z^2/2 - eps^2/(2 sigma_v^2) as
the difference of two terms that each reach 1e14 once sigma_v nears its
floor. The rounding noise came out as a spurious gain, which made the
degenerate sigma_v -> 0 fits look better than they were. That difference is
now formed exactly (for NG it equals eps/sigma_u + sigma_v^2/(2 sigma_u^2),
for NNAK -mu eps^2 / sigma^2), with the remaining factor supplied by a new
.log_pcf_w() that is exact in the far tail. Both likelihoods now equal the
convolution integral to 1e-8 at sigma_v down to 1e-7. The efficiency
predictions, which formed the same terms, are changed the same way.
Third, with the likelihood exact, NNAK's surface often has a higher mode at
small shape that its single half-normal moment start does not reach. NNAK
now goes through the same candidate-and-polish multistart as NG, with the
shape swept down to 0.01. A second family of candidates covers a mode
near shape 0.001 in which u is mostly near zero with a long tail: for
shapes 0.001 to 0.03 it sets sigma_u^2 to 10, 20 or 40 percent of the
residual variance rather than anchoring E[u] on the residual skewness,
which puts sigma_u near 0.03 there. The two families are ranked apart,
so the second adds to what is polished rather than displacing the first.
Over 216 fits (12 samples by 18 response columns), the number of NG fits ending more than one log-likelihood unit below the NE it nests falls from 7 to 0 (worst shortfall 217.4 to 0.02), and the number of NNAK fits ending that far below NHN falls from 18 to 0 (worst 15905.5 to 0). No NNAK fit ends lower than on main by more than 0.014, and 93 end higher. NNAK fits take about 65 percent longer.
sfm(model_name = "TSL") could return a non-physical fit. Both of the
truncated skew-Laplace likelihood's exponentially tilted Gaussians were
formed as two numbers whose large parts cancel analytically: with
t = sigma_v/sigma_u the tilt raises t^2/2 + w t and log Phi takes
exactly that back, so A - aa^2/2 = -eps^2/(2 sigma_v^2). As sigma_u
collapsed, t grew without bound and the two pieces cancelled
catastrophically -- while both stayed FINITE, so the non-finite guard never
saw it. At sigma_v = 1.126e20, sigma_u = 1e-7 the old form returned a
log-likelihood of +2.55e40, roughly 1e40 better than any legitimate
optimum and therefore an attractor. Both tilts, and the matching pair in the
TSL efficiency block, now go through .log_phi_tilt(), as NE and NGE
already did.
Over 432 TSL fits spanning 24 samples and 18 response columns -- mostly
deliberately misspecified, which is where sigma_u collapses -- 8 returned
a non-physical log-likelihood, and 8 more returned a perfectly plausible one
that was wrong by 1 to 11 log-likelihood units. Both groups are now correct.
The second group is the one worth knowing about: those fits reported values
like -353.45 where the right answer was -364.47, with nothing about the
number itself to suggest a problem.
Every one of the 16 sat at sigma_u = 1e-7 exactly -- its floor -- while
scoring ABOVE the OLS that TSL nests at sigma_u = 0. That conjunction is a
usable diagnostic: it flags 17 of the 432 fits before the fix and 0 after.
Reported as issue #28.
sfm(model_name = "tHN") lost every standard error on about a third of
fits. The tHN branch returned -.Machine$double.xmax / n per observation
for a parameter outside its domain. Summed over the observations that is
.Machine$double.xmax, and optim() differences the objective to build its
gradient, so the first difference was not merely large but non-finite, and
the Hessian -- and with it every standard error -- came back NA.
The guard looked unreachable: nu is bounded below at 2.05 while the guard
tests nu <= 2, and sigma_u at 1e-7 while the guard tests
sigma_u <= 0. It is reached because optim() steps outside its own box
to form the finite difference -- ndeps defaults to 1e-3, and the guard
was observed firing at sigma_u = -0.0009999, exactly 1e-7 - 1e-3. Over
nine fits, every fit that reached the guard lost all its standard errors and
every fit that did not kept them.
Over 84 tHN fits the number losing every standard error falls from 25 to
2, with the fitted log-likelihoods unchanged. The two that remain sit at
nu = 7720 and nu = 69 -- the t's normal limit, where the degrees of
freedom are weakly identified and the Hessian is legitimately near-singular.
The penalty is now .SFA_CONSTANTS$DOMAIN_PENALTY_PER_OBS, chosen by
measurement rather than by analogy: the obvious fix of matching the 1e12
used by the NGE, NGB2, NLN and NW branches recovers nothing, since
1e12 differences just as badly. Those four branches are unchanged here --
whether their own guards are reachable the way tHN's is has not been measured,
and assuming it would repeat the reasoning that left tHN untouched.
A covariance that failed on a badly scaled Hessian now succeeds, and one
that genuinely cannot be computed now names the parameter responsible.
Two separate defects sat behind the same "no invertible Hessian" message.
First, solve() was called on the raw Hessian, whose diagonal can span
twenty orders of magnitude when one scale parameter is pinned near zero and
the likelihood is correspondingly sharp in that direction; the matrix is
then declared computationally singular even though it is well conditioned
once equilibrated. vcov() and the sandwich methods now solve
D^(-1/2) H D^(-1/2) and undo the scaling, which is algebraically the same
inverse and numerically a far better one -- on the fit that prompted this,
rcond rose from 3e-17 to 8e-2 and the standard errors appeared. Second,
when no inverse exists on any scaling, the error said only that one did not,
leaving the user to guess which of eleven parameters had collapsed. It now
reports the specific cause per parameter: a score column that is
identically zero, a coefficient pinned on its bound, a flat direction of the
diagonal-normalised Hessian, or non-positive curvature.
psfm_bootstrap()'s seed_offset did not give distinct runs. Replication
b was seeded with b + seed_offset, so a run at seed_offset = 1 repeated
replications 2 to BOOT of the run at seed_offset = 0 exactly, although
the help page promised "reproducible-but-distinct seeds across multiple
bootstrap runs". Two runs pooled as independent shared all but one draw, and
a standard error from the pool understated the spread. The seed is now
b + seed_offset * 100000, so different offsets never share a seed for
BOOT up to 100000. seed_offset = 0, the default, seeds exactly as
before, so existing results reproduce; results from a non-zero offset
change. Measured alongside: the draws are identical whatever numCores is
and whatever the caller's RNG kind, and the caller's random stream is left
as it was -- both now pinned by tests.
sfm("NG") evaluated its log-likelihood as exactly 0 near a collapsed
scale. The density contains z^2/4 + log D(z), the
parabolic cylinder function's log plus a quadratic, and the log-D routine
returns -z^2/4 plus a term of order log z. Added as two numbers they
cancel, and for NG z is about sigma_v / sigma_u, so the cancellation
bites exactly where sigma_u approaches its lower bound. On 200
observations the summed NG log-likelihood was off by 0.03 at
sigma_u = 1e-7, by +117 at 1e-9 -- spuriously better than the truth --
and exactly 0 from 1e-10 down. Zero is finite, so the non-finite guard let
it through, and it sat far above any real optimum. The multistart guard
did not always catch it: a plain sfm(y ~ x1 + x2, model_name = "NG") on
data_gen_cs(N = 300, rand = 12, sig_u = 0.5, ...) returned an objective
of exactly 0 at sigma_u = mu = 1e-7, with no warning, and now returns an
interior fit whose log-likelihood beats the nested normal regression's. The sum is now formed as one
quantity, by the large-argument expansion where it would cancel, and the
likelihood tends to its correct limit, the normal density of v alone.
"NNAK" contains the same combination and now forms it the same way, but
was measured not to be exposed: its argument grows only as sigma_v -> 0,
where the likelihood is already, correctly, very negative.
The same combination made efficiency()'s NG and NNAK scores NaN once
z passed about 53, as exp(z^2/4) overflowed in numerator and
denominator alike. That ratio is now formed in logs.
optHessian = FALSE with the default PSopt = FALSE failed outright for
sfm()'s default model and for every lcsfm() model, and lost the
log-likelihood everywhere else. Without the final optim() stage a fit
stored the earlier stage's result as $opt unchanged, and two of those
stages do not use optim()'s field names: nlminb() reports objective,
bobyqa() reports fval and ierr. Every consumer reads $value.
sfm()'s "NHN", "NE", "NTN" and "NU" run nlminb instead of bobyqa,
and nlminb's result was discarded, so $opt was NULL and the fit stopped
with "number of items to replace is not a multiple of replacement length".
The other eleven sfm() models returned a fit whose logLik(), AIC and BIC
were NA, with a warning that blamed moment-based estimation. lcsfm()'s
"LCM" and "LCM_Z" stopped with "invalid argument to unary operator",
and "LCM_CN" never kept any stage's result at all. zsfm(), ttsfm(),
selsfm(), copsfm(), ivsfm() and psfm() had the silent-NA half.
An earlier stage's result is now put into optim()'s shape before it
becomes $opt, and sfm() keeps nlminb's result when that is the stage
that ran. logLik() on such a fit is the log-likelihood at the reported
estimate. optHessian = TRUE (the default) and PSopt = TRUE are
unchanged.
fitted(), residuals(), predict() and efficiency(logDepVar = FALSE)
on a zsfm() or ttsfm() fit could still answer from the wrong data. The
1.2.1 check on recovered data (A47) is decisive where a fit stored its OLS
residuals, as sfm() does: the recovered frame is refitted and compared.
zsfm() and ttsfm() stored none, so for them the check could only count
rows -- and ttsfm("TTNLS") stored no row count either. A different data
frame of the same shape, bound to the name the call recorded, passed, and
every method answered from it without a warning. Both entry points now keep
the frontier's OLS residuals on the fit, and a same-shaped decoy is refused.
ttsfm("TTNLS") also now stores nobs.
ttsfm("TTHN") replaced optim()'s result on a tie, and warned that the
final stage had failed. The branch keeps whichever optimizer stage
attained the lowest objective, and chose with a bare which.min(). On one
seed bobyqa and optim() had converged to the same point -- parameters
agreeing to 1e-10, identical Hessian eigenvalues -- and a difference of
2.8e-14 at an objective of 189, one unit in the last place, discarded
optim()'s result. The warning printed both objectives to four decimals,
where they were equal. optim() is now displaced only by a stage that
improves on it by more than optim()'s own convergence tolerance
(reltol = sqrt(.Machine$double.eps), relative), and the warning prints
the objectives to ten significant figures. Reported estimates change only
where the old choice was a tie, and then by the width of the tie. Standard
errors can change more on such fits, because keeping optim() also keeps
optim()'s Hessian, where the tie used to hand over a numDeriv one. On
one seed the SE count went from 6 to 4. That fit had collapsed sigma_v
to 0.001 (the corner described in ?ttsfm), where the objective is flat
to about 5e-4 and no finite-difference Hessian resolves the curvature, so
neither count was meaningful there. What changes in general is that the
SEs no longer depend on which stage won a tie.
sandwich::sandwich() and sandwich::vcovCL() returned silently wrong
standard errors for every model whose estimation and reported parameter
scales differ. bread() was built from vcov(), on the REPORTED scale,
while estfun() returns scores on the ESTIMATION scale; sandwich composed
the two without knowing they disagreed. Where the two scales have the same
dimension this produced finite, plausible, wrong numbers: on
ttsfm("TTNE") the three sigma standard errors came out 29% low, 18% low
and 26% HIGH -- each wrong in a different direction -- while the betas,
whose Jacobian is 1, were correct. Where the dimensions differ it could not
conform at all, which is how ivsfm("IVLIML") (11 estimated, 6 reported)
surfaced it, failing with non-conformable arguments.
Affected ttsfm(), copsfm(), ivsfm() and psfm()'s PL80 family.
sfm(), zsfm(), lcsfm() and psfm()'s TRE/GTRE family report what
they estimate and were never affected.
bread() is now on the estimation scale, matching estfun(), so
sandwich::sandwich() and sandwich::vcovCL() return an internally
consistent estimation-scale covariance -- and are defined for IVLIML at
all. For a covariance on the scale coef() reports, use the new
vcov(type = "sandwich") and vcov(type = "clustered"), which apply the
mapping. Where the two scales agree the results are unchanged.
psfm(model_name = "GTRE") combined its antithetic simulation draws on the
log scale. The persistent component is simulated as +r and -r, two
draws of the same density, which average as log((A + B) / 2). Both call
sites -- the loop path and the vectorised one -- computed
0.5 * (log A + log B), which is log(sqrt(A * B)): the geometric mean,
which AM-GM puts strictly below the arithmetic one whenever the halves
differ. The vectorised site's own comment said the density was being
averaged.
This is a finite-draw bias rather than an asymptotic error -- each half is
separately a consistent simulator, so the two rules agree as the draw count
grows. But the default count is ceiling(sqrt(nrow(data))) + 100, about 118
draws for a 60x5 panel, and measured there on identical draws the summed
log-likelihood was 0.23 too low, with the sigma_r profile optimum shifted
down 4.8%. The bias falls to 0.012 at 3200 draws and 0.0013 at 20000.
Averaging is now done on the density scale through log-sum-exp, so neither
half is exponentiated on its own.
GTRE is the only model affected: TRE uses the +r half alone and
GTRE_Z has no antithetic pairing. GTRE estimates move slightly as a
result. Raised by Chris Parmeter.
The optimizer's credibility gate accepted points outside the feasible
region. .sfa_point_credible() decides whether an earlier optimizer stage
reached a point worth re-polishing from. It correctly refused a point resting
ON a parameter bound -- the signature of a scale driven to its floor -- but
tested that with abs(par - lo) <= tol, which measures distance rather than
direction. A point sitting BELOW the floor is far from it by that measure and
so was admitted, meaning the gate refused the boundary and accepted the region
beyond it. psfm("TRE_Z") is how this surfaced: sigma_r = -0.2 against a
floor of MIN_POSITIVE reads as 0.2 away from the bound, so a negative
standard deviation passed as a credible optimum. Points outside their bounds
in either direction are now rejected.
psfm(model_name = "TRE_Z") could report a NEGATIVE sigma_r. The
sign of sigma_r is not identified: it enters the likelihood only as a
scale on simulation draws that are symmetric about zero, so +sigma_r and
-sigma_r describe the same model and fit equally well (flipping it moves
the objective by 8e-6 in relative terms; flipping sigma_v, whose sign is
pinned by sigma_u / sigma_v, moves it by a factor of 7). TRE_Z was the
only psfm() model that gave bobyqa an unbounded lower limit instead of
a per-parameter vector flooring the two scale parameters, so it was the
only one that could settle on the negative root -- about 7% of samples,
against none for TRE, GTRE, GTRE_Z or GTRE_FML. Every later
optimizer stage floors sigma_r at MIN_POSITIVE, so once bobyqa
returned a negative value it was already outside their box and could not be
corrected. sigma_r is now floored from the first stage. Affected fits
return the positive value of the same magnitude; a sample with no
detectable random effect goes to the boundary rather than to a negative
number.
This also removes a spurious convergence failure. The convergence/
harness read sigma_r's MSE as flat in the sample size -- slope +0.010 at
R2 0.008, where every other parameter of the same model passes at root-n --
because roughly one fit in fourteen contributed the squared distance
between +sigma_r and -sigma_r regardless of n. On absolute values the
same archive gives slope -0.949. The estimator was consistent all along.
estfun() could difference against a numerical-stability penalty and
return scores of order 1e150. When a finite-difference step leaves a
likelihood's admissible region the closure returns a large penalty as a
sentinel. That penalty is finite, so it passed the existing
is.finite() guard and was differenced like an ordinary value, filling the
meat matrix with enormous numbers that looked like numbers. Any fit resting
at or near a parameter bound was affected, which is the worst case for it to
reach: the OPG is offered precisely because it needs no Hessian and so works
where the default path fails. Steps that leave the region now fall back to a
one-sided difference on the admissible side, are zero only where both sides
are inadmissible, and raise a warning saying the fit sits at a bound and
that OPG/BHHH errors understate uncertainty in those directions.
sfa_diagnostics() labelled estimation-scale quantities with reported
names. The Hessian, the numerical gradient and plot()'s likelihood
slices are all indexed by opt$par, which for ivsfm("IVLIML") has 11
entries against 6 reported names. The gradient came back as
(Intercept), x1, x2, sigma_u, sigma_v, rho_x2, NA, NA, NA, NA, NA --
entries 4 to 6 carrying names that belong to different quantities, which is
worse than no name at all. These now use the estimation-scale names, and
print() discloses the second scale when it differs. The correlation
diagnostic deliberately keeps the reported names, because it is built from
vcov(), which is already on that scale.
sfa_diagnostics()$convergence$logLik went through -opt$value and so
reported a divergence for robust fits; it now goes through logLik(). The
achieved objective is still reported, as objective, under a name that
says which objective it is.
estfun() gave robust fits advice that could not help them. The
"no retained likelihood" check ran before the robust check, so a robust fit
made without keep_objective = TRUE was told to refit with it -- and the
refit then hit the robust refusal. The robust check now runs first.
logLik(), AIC() and BIC() returned the divergence objective for
robust fits. sfm(robust = "mlqe" | "psi" | "mdpd") minimises an MLq,
psi or density-power divergence, not a negative log-likelihood, so
-opt$value is not a log-likelihood and criteria built on it are not
comparable with anything. On a 300-observation NHN fit the numbers were:
| fit | logLik() | AIC() |
|---|---|---|
| mle | -299.50 | 608.99 |
| mlqe | -256.77 | 523.53 |
| psi | +29408.72 | -58807.43 |
| mdpd | +1240.59 | -2471.18 |
Comparing AIC(mle) against AIC(robust) -- the obvious thing to do --
made the robust fit look better by 3080 points, entirely from the change of
objective, with no warning. mlqe is the most dangerous of the three
because -256.77 sits plausibly beside the MLE's -299.50. estfun() and
TIC() already refused these fits in explicit terms; logLik() now
refuses in the same terms, returning NA with a warning and keeping its
df/nobs attributes so AIC()/BIC() propagate NA rather than error.
sfma() was never exposed: it passes robust to every candidate, and
since robust supports "NHN" only, the others fail and it stops first.
print() and summary() labelled that same number "log likelihood".
They now name the objective they are showing -- mdpd objective (NOT a log likelihood): for a robust fit, penalised log likelihood: when
penalty_c > 0 -- and are unchanged for ordinary maximum-likelihood fits.
logLik() reported the penalised objective for lcsfm() fits with
penalty_c > 0, and the field meant to prevent that was itself wrong.
like.fn() returns -(L + P) for log-likelihood L and penalty P, so
-opt$value is L + P; the plain log-likelihood is -opt$value - P. Both
the LCM_CN and the LCM/LCM_Z branches stored -opt$value + P, giving
L + 2P -- so the "unpenalised" figure was further from the truth than
the penalised objective it was meant to correct. P is at most zero (by
Jensen, J log J + sum log pi <= 0), so both understated the
log-likelihood and AIC()/BIC() over-penalised the fit.
Measured on a 600-observation two-class mixture with 74/26 shares at
penalty_c = 20: logLik() returned -804.49 where the log-likelihood is
-799.35, understating by 5.14 and overstating AIC() by 10.28. The error
is zero at penalty_c = 0 and grows with the penalty and with class
imbalance, so it is invisible in the balanced case.
logLik() now reads the corrected field, and the LCM_CN penalty is
evaluated at the parameters the fit actually reports. Verified against the
definition: logLik() equals sum(objective(par, per_obs = TRUE)), the
sum of the unpenalised per-observation contributions, at every
penalty_c tested.
lcsfm_homogeneity()'s MLR statistic is unaffected -- it uses the
penalised objective deliberately, since the penalty is what makes the
likelihood ratio "modified". Only the reported diagnostic changes.
AIC(), BIC() and TIC() under-penalised ivsfm("IVLIML"). Its
degrees of freedom came from length(coefficients), the number of
parameters reported. IVLIML maximises over reduced-form nuisance
parameters it never reports -- 11 estimated, 6 reported -- so it was charged
for 6, a free 2 * 5 = 10 AIC points against any model that reports what it
estimates. On a 300-observation sample whose two log-likelihoods differ by
0.005, AIC() preferred IVLIML by 0.010; with the correct df it prefers
IVCF by 9.99, so the defect was large enough to reverse the choice.
TIC(detail = TRUE) inherited it in a second place. Its penalty is a
trace computed on the estimation scale, so comparing it against the
reported count gave ratio = 11.898 / 6 = 1.98, which reads as a gross
violation of the information matrix equality. Against the estimated count
it is 11.898 / 11 = 1.08 -- what a correctly specified model should give.
df is now the estimated count wherever the fit carries an optimizer,
falling back to the reported count for the moment-based and FE estimators
that do not. Every other model in the package reports exactly what it
estimates, so nothing else changes.
vcov() was permuted for psfm("PL80"), "BC92", "K1990" and
"K1990modified". These models estimate (sigma_v, sigma_u, beta) and
report (beta, sigmaSq, gamma), so the two scales differ by a permutation
as well as a transformation. vcov() returned the estimation-scale inverse
Hessian under the reported names, which put Var(sigma_v) on the intercept
and Var(beta) on sigmaSq. On a 60-firm panel confint() reported the
intercept as [1.0101, 1.0707] where the printed standard error implies
[0.9303, 1.1505] -- 3.6 times too narrow, with every other row wrong in a
different direction. The delta-method Jacobian these fits already computed
for their own standard errors is now stored on the object and applied by
vcov(), so the two agree exactly.
psfm("PL80_MVTN") stored out untransposed. Every other entry point
stores the p x 3 transpose, and the documented access fit$out[, "par"]
failed on this model alone with subscript out of bounds.
psfm("PL80_MVTN") could fail to start, and could stop short of a
maximum. Its out-of-domain penalty was .Machine$double.xmax; it is the
one likelihood handed straight to optim() with no derivative-free stage
in front of it, so L-BFGS-B differenced that penalty against an ordinary
value and overflowed, aborting the fit with non-finite value supplied by optim on some samples. Separately, the orthant probability inside the
likelihood is computed to limited precision, so the finite-difference
gradient is noisy and a single L-BFGS-B run stopped at points that were not
maxima while still reporting convergence = 0. The penalty is now damped
and a derivative-free stage runs between two L-BFGS-B runs, as in every
other psfm() model. Already-converged fits are unchanged to the last
digit; where the single stage quit early the fit gains 5 to 7
log-likelihood points.
vcov() returned standard errors on the wrong scale for ttsfm(),
copsfm() and ivsfm(), and failed outright for ivsfm("IVLIML"). These
entry points estimate log sigmas and report sigmas, and ivsfm() reports
rho = t/sqrt(1 + t't). Their reported standard errors have always carried
the delta-method correction, but vcov() did not: it named the
estimation-scale inverse Hessian with the reported names and returned it. So
sqrt(diag(vcov(f))) disagreed with the sigma_u standard error the same
fit printed by a factor of exactly 1/sigma_u -- wrong in a way that looks
entirely plausible, and inherited by confint(), lmtest::coeftest() and
anything else built on vcov(). ivsfm("IVLIML") additionally estimates
reduced-form parameters it never reports, so naming an 11x11 matrix with 6
names failed with length of 'dimnames' [1] not equal to array extent.
vcov() now carries the fit's own Jacobian from the estimation scale onto
the reported one, and marginalises estimated-but-unreported parameters by
taking the corresponding block of the inverse rather than inverting the
block -- the latter conditions on the nuisance parameters being known and
understates every standard error. Both the Hessian and the BHHH path go
through the same map, so the two can no longer disagree about scale.
influence_sfa() was pairing estimation-scale scores with a reported-scale
covariance for the same reason, and is fixed with it.
ttsfm(model_name = "TTHN") could stop with 'x' must be an array of at least two dimensions. When the three optimizer stages disagreed about
which reached the best objective, ttsfm() substituted the winning stage's
result into the fit -- but a bobyqa() or psoptim() object does not carry
the hessian that the standard-error step then asked for, so the fit died
after the work was done. The substituted result is now given optim()'s
shape, with a Hessian computed at its own parameters. The same substitution
also silently cost the fit its log-likelihood: bobyqa() reports the
objective as fval where optim() reports value, so a bobyqa-won TTHN
fit returned NA from logLik(), AIC() and BIC() with no warning.
Both are fixed together.
One consequence worth knowing. Where the stages tie -- and on this
likelihood they frequently tie to within rounding -- which stage "wins" is
decided by the last bits of the objective and can differ between platforms.
The substituted path computes its Hessian with numDeriv, while an
unsubstituted fit carries optim()'s own. On an ill-conditioned fit the two
estimators can disagree about whether the Hessian is positive definite, and
therefore about how many standard errors are reportable. This is a property
of the conditioning, not of the fix -- but before the fix these fits errored
rather than reporting anything at all.
Cost frontiers could be read as production frontiers by three
post-estimation functions. composed_cdf(), lcsfm_homogeneity() and the
robust residual accessor each decided orientation their own way. Two were
wrong in the same manner: a fit made with inefdec passed as a variable
rather than a literal -- flag <- FALSE; sfm(..., inefdec = flag) -- had the
unevaluated symbol forwarded or an unresolvable value silently treated as a
production frontier. A cost fit was then scored with the production sign, with
no error and no warning. All three now resolve orientation through the same
helper the rest of the package uses, and stop with a message when the
orientation genuinely cannot be determined rather than assuming one.
efficiency(logDepVar = FALSE) gave cost frontiers the production
formula. On the level scale it returned 1 - u/f for every fit; for a cost
frontier, where y = f + u + v, efficiency is minimum over actual cost,
f / (f + u), the level-scale counterpart of the exp(-u) the default log
scale already reports. Cost-frontier level-scale scores from 1.2.0 were too
low (negative where u exceeded f). Production fits are unchanged. The
orientation is now read from the fit's inefdec even when it was passed as
a variable, which predict(type = "response") also could not do. Reported
in GitHub issue #2.
efficiency_ci() could pair point estimates with the wrong
observations if a fit's stored u_hat or exp_u_hat had a length that
divides the number of observations: a data-frame column silently recycles
such a vector. It now stops with a message instead (issue #2).
sfm(model_name = "NG") and "NNAK" could return a fit far worse than a
point they had already reached. The three optimizer stages start from the
best of NG's polished multistart candidates, and could end below it: on one
simulated sample the start scored -259.6 in log-likelihood and the fit came
back at -449.5. The fit is now never worse than that start; where the stages
end below it, the final stage is re-run from it and the better result kept.
Over 40 simulated samples this changed two fits by about 180 log-likelihood
points and left the rest within 0.5.
density_weights() gave the same weight to every observation whose fitted
density fell below the smallest positive double, and lost the ordering among
them. The weight was formed by evaluating the density on the natural scale,
where .dens_nhn() floors at .Machine$double.xmin, and then raising it to
the power c. Every observation past that floor -- log f below -708.4 --
arrived at the floor itself and came back with one identical weight. The
effect is largest exactly where the diagnostic is used: on a sharply peaked
fit (sigma_v near 0.04, lambda near 24, as maximum likelihood returns on
the coal data of Bernstein, Parmeter and Wright (2026)), 63 of 3511
observations shared a single weight of 1.7e-67 while their true log densities
spanned -4680 to -708, and order(w) -- the documented way to ask which
observations an estimator discounts most -- was wrong among them. The weight
is now formed as exp(c * log f), normalised in logs before exponentiating.
A new log argument returns log w directly, which is what ranking or
thresholding should use: the smallest weights here run to exp(-1042) and
are not representable, let alone distinguishable, once exponentiated.
Weights that were already above the floor are unchanged to about 1e-14, and
no fitted value moves: hscore() and the robust objectives are bit-identical
before and after, and calibrate_c() shifts by about 1e-12 relative.
The robust objectives formed exp(loglik)^c rather than exp(c * loglik).
.robust_objective_vec() and the calibration's .per_obs_term() exponentiated
the log-likelihood before applying the power, so an observation with loglik
below about -745 contributed exactly the c-limit of its term instead of a
value near it. This is now formed in the exponent throughout. No estimate
changes -- at those depths the recovered term is below the resolution of the
1 it is subtracted from, so the contribution is still 1/c to the last bit --
but the objective is no longer discontinuous at the underflow boundary, which
matters for the numerical derivatives taken through it.
psfm(model_name = "FD"), "TFE" and "TFE_WMLE" were not reproducible,
and changed the caller's random-number stream. Their variance components are
seeded by an internal pooled fit whose own starting values were unseeded
runif() draws. Two fits of identical data in one session could differ by
up to about 3e-4 in relative terms, and occasionally one of them stopped with
an abnormal line-search termination where the other did not; every fit also
advanced the user's RNG. The draws are now taken under a fixed local seed and
the caller's RNG state is restored. The starting-value distribution is
unchanged, but results will differ slightly from fits made with 1.2.0, which
were themselves random.
psfm(model_name = "GTRE") by simulated ML, and "GTRE_Z", reported
efficiency scores that were not reproducible, and changed the caller's
random-number stream. The persistent and transient efficiencies (H,
U) are ratios of multivariate normal orthant probabilities from
tmvtnorm::ptmvnorm(), which integrates by randomized quasi-Monte Carlo.
Two fits of identical data differed by up to 0.4% in U and 0.2% in H.
Parameter estimates and standard errors were never affected, nor were
estimator = "fiml" and "TRE". The integration now runs under a fixed
local seed and the caller's RNG state is restored; the scores still carry
the integrator's own error of about 1e-3, now the same error every time.
ttsfm(model_name = "TTNE") reported a wrong Eemu.cond and M6 metric
whenever sigma_v was not 1. The exponent in E[exp(-u) | eps] used
0.5 * sig.v where the closed form has 0.5 * sig.v^2. At sigma_v = 0.3
the reported value was 10.9% too high (0.522 against an exact 0.471). The
square is restored; all six conditional expectations now agree with numerical
integration to about 1e-11. Parameter estimates were not affected.
ttsfm(model_name = "TTNE") clipped its log-likelihood's exponents. The
density was formed as log(pnorm(beta) * exp(alpha) + pnorm(b) * exp(a)) in
levels, with alpha and a capped before exponentiating to avoid overflow,
so the reported density was wrong wherever the cap bound and either pnorm()
could underflow. It is now a log-sum-exp of the two terms, exact with no cap.
On 20 simulated two-tier samples, 19 fits are unchanged to about 1e-6 and one
reaches a log-likelihood 1.2e-4 higher.
ivsfm() could report efficiencies above 1. The JLMS predictor divided
dnorm(zz) by pnorm(zz) floored at machine epsilon, so for a firm well above
the frontier (zz below about -8.1) the inefficiency estimate came back
negative and exp(-jlms) exceeded 1: on simulated data with one firm 20 above
the frontier, model_name = "C2SLS" reported an efficiency of 71, and 25
million at 80 above. It is now computed in logs and cannot be negative.
Estimates are unaffected.
zsfm(model_name = "ZISF_Z") often stopped in a spurious local optimum.
Its likelihood has a second basin in the regime-link coefficients, 10 to 64
log-likelihood units below the best maximum. In 50 simulated samples (40%
of firms efficient, the covariate unrelated to the regime) 7 shipped fits
ended there, reporting a regime slope averaging 0.38 against a truth of 0,
with no warning. "ZISF" already guarded against the same problem with
several starting fits; "ZISF_Z" now does too, over the link intercept. On
the same 50 samples no fit falls more than 0.53 below the best of six
polished starts, none is worse than before, and a fit takes about 0.20 s
instead of 0.13 s.
zsfm() and lcsfm() could return NaN or -Inf for an observation far
from the frontier. zsfm()'s inefficient-regime log-density was
log(dnorm()) + log(pnorm()), which is -Inf once either factor underflows
(beyond about 38 standard deviations; the truth there is finite, e.g. -919),
and both zsfm() and lcsfm() computed the JLMS predictor as
dnorm(zz) / pnorm(zz), which is 0 / 0 below zz of about -37 (a firm far
above the frontier; the truth is about 0.007). Both are now taken in logs. On
ordinary data the likelihood is unchanged at every fit checked, ZISF and
LCM fits agree to within 7e-5, and lcsfm() is identical.
sfm(model_name = "NU")'s log-likelihood, and "NHN"/"NHN_Z"'s at
extreme parameters, were floored rather than computed. NU took
log(pnorm(b) - pnorm(a)) as a difference floored at .Machine$double.xmin:
for a residual more than about 9 noise standard deviations above the frontier
both CDFs round to 1, and the log-density came back as -708 where the truth
is -44. NHN and NHN_Z floored the density in levels, so any log-density below
-708 was reported as -708. Both now work in logs. Estimates move by less than
8e-8 (NHN), 5e-6 (NHN_Z) and, in flat cases, 3e-3 (NU) at the same maximum
likelihood.
Battese-Coelli efficiencies could come back as 0 for firms well above the
frontier. sfm()'s "NHN" and "NHN_Z" and psfm()'s "TRE" and
"TRE_Z" computed exp_u_hat (U) as a ratio of 1 - pnorm() terms, which
underflow together once the scaled residual lambda * eps / sigma passes
about 8: at lambda = 10 a firm 0.8 above the frontier was reported at 0.83
instead of 0.988, and one 1.0 above at 0 instead of 0.990. A single NHN fit
tends to keep this out of reach, because such a point pulls lambda down; the
exposed cases are "NHN_Z" and "TRE_Z", where lambda varies by observation.
They now use the log-space predictor the package already used elsewhere;
ordinary fits agree with the old values to 3e-15. meanefficiency()'s
closed forms for half-normal, truncated-normal and Rayleigh inefficiency,
which overflowed to NaN for a scale or truncation point near 40, are now
evaluated in logs too.
sfm(model_name = "NR") was drawn to sigma_v = 0 by a rounding artefact.
Its log-density formed -(eps/sigma_v)^2/2 + z^2/2, two very large terms
whose exact sum is -eps^2/sigma^2. As sigma_v fell toward zero the
difference lost its digits and read too high (by 22 log-likelihood units per
observation at sigma_v = 1e-9), so fits ran to the sigma_v bound and could
report a log-likelihood far above the true one (-539.9 against -585.1 on one
simulated sample). The term is now formed exactly. On data with
lambda = 2 nothing changes; with lambda = 10 fits reach the best maximum
of the exact likelihood in 187 of 200 samples.
TIC() could return an absurd value for a fit on a parameter bound, and
sfma(weights = "tic") then gave that model all the weight: -2.64e14 for an
NR fit with sigma_v on its bound, although another candidate fitted
better. TIC() now stops when the Hessian is not positive definite or the
Takeuchi penalty is not positive; sfma() gives such a candidate no weight
and says which one and why (unless quiet = TRUE), and vuong() refuses it.
sfm(model_name = "NR") evaluated its log-likelihood inaccurately for
residuals far above the frontier. The term log(phi(z) - z * Phi(-z)) was
computed as the difference of two nearly equal quantities through erf():
from z of about 7 it lost all precision and overstated the log-density by up
to about 7 (by 0.8 at z = 8 and 4.6 at z = 10). Such z values arise at trial
parameters with a small sigma_v, so the optimizer was working on a
distorted surface there, and exp_u_hat used the same form. Both are now
computed through the continued fraction for the Mills ratio, accurate to
1e-13 across the range. On a design with lambda = 2 every fit reaches the same
maximum as before; with lambda = 10 fits reach the best maximum as often as
before (156 against 155 of 200) but can differ individually where the
likelihood is flat.
psfm(keep_objective = TRUE) was silently ignored for every model except
"GTRE" (simulated ML), "TRE", "GTRE_Z" and "TRE_Z", and estfun()'s
error then advised refitting with keep_objective = TRUE -- which could not
help. psfm() now warns when the argument has no effect, and the message
says which models retain their likelihood.
psfm(model_name = "GTRE_Z") lost its whole fit when one firm's efficiency
could not be computed. A single firm whose posterior system could not be
inverted stopped the fit with "Failed to compute GTRE posterior matrices",
discarding estimates that had already converged. That firm's U and H are
now NA, with a warning naming it.
psfm(model_name = "GTRE_Z") stopped on some platforms with "sigma must be
positive definite" from tmvtnorm, when a firm's efficiency posterior
covariance lost positive definiteness to rounding: the same call fitted on
macOS and failed on Windows and Linux. The covariance is now repaired when
the defect is rounding-sized (eigenvalues clipped at a relative floor, for
"GTRE" by simulated ML as well); a larger defect gives that firm NA
efficiencies with the warning above. Where the covariance was already valid
the results are identical.
A starting value just outside its bound stopped the fit with "Starting
values violate bounds" from bobyqa(). On a panel whose random-effects start
puts the firm-effect scale at essentially zero, psfm(model_name = "TRE")
began with sigma_r a hair above its 1e-7 bound on one platform and below
it on another, so the same call fitted on macOS and failed on Linux and
Windows. The internal bobyqa() stage now moves only coordinates that lie
strictly outside their bounds to just inside; every start already inside is
used unchanged, so fits that ran before are identical.
psfm(model_name = "GTRE_Z") and "TRE_Z" stopped with
"0 < ctrl$rhoend is not TRUE" when the random-effects regression that seeds
them put the firm-effect variance on its boundary. Its firm effects are then
all zero, and the intercept-only fit that turns them into starting values was
started at zero scale, which leaves bobyqa() no trust region. It now starts
from a small positive scale; inputs with any spread are handled exactly as
before.
psfm() died with "system is exactly singular" for PL80, BC92,
K1990, K1990modified and SSFE on a between-rank-deficient design,
such as factor(year) in an unbalanced panel. None of these models uses the
random-effects starting regression -- each builds its own starting values --
but start_panel() ran it for them anyway, and did so outside the
collinearity guard, which skips these models. The regression is now run only
for the seven models that read its result (TRE_Z, GTRE_Z, TRE, GTRE,
GTRE_FML, GTRE_SEQ1, GTRE_SEQ2). Reported from a replication of
Filippini and Greene (2016) on the Swiss railways panel, which had been
working around it by substituting a different variance-component method
inside plm().
psfm(model_name = "SSCRE") could not fit an unbalanced panel, and
most real panels are unbalanced. In the between-firm dimension each Mundlak
mean duplicates its regressor, so plm's Swamy-Arora variance components
for an unbalanced panel hit an exactly singular system. SSCRE's variance
components are now defined as the Swamy-Arora components of the model
without the Mundlak means. On a balanced panel that is exactly what it
computed before, so balanced results are unchanged; on an unbalanced one the
identical components come from plm's Wallace-Hussain step on the augmented
model, checked to agree. The slopes on time-varying regressors still equal
the within estimates exactly.
psfm(model_name = "SSRE") and "SSCRE" now explain a singular
random-effects fit instead of failing with a bare "Lapack routine dgesv:
system is exactly singular". When the formula is collinear between firms the
message names the columns and suggests remedies.
pracma and MASS are no longer imported, and pbapply moves from
Imports to Suggests. Each was used for a single function. erfinv() and
ginv() are now internal copies that compute the same thing the same way:
the simulated-ML draws and every fit are bitwise identical. pbapply only
ever supplied psfm_bootstrap()'s optional progress bar, and the code already
ran without it; its pkgs default no longer loads it on the workers.
tmvtnorm is replaced by mvtnorm, which tmvtnorm itself was built on:
13 recursive dependencies become none outside base R. psfm()'s GTRE
efficiency scores called tmvtnorm::ptmvnorm(), which for an untruncated
orthant is exactly mvtnorm::pmvnorm(); they now call that directly and are
bitwise identical. data_gen_p()'s truncated multivariate normal draw for the
PL80_MVTN column now comes from an internal Gibbs sampler (checked against
the exact truncated moments), so that column and y_pl_mvtn differ from
earlier versions for the same seed; every other column is unchanged.
gsl is no longer imported (it moves to Suggests, where the tests use
it as an independent reference), so installing sfa from source no longer
needs the GNU Scientific Library. Its only use was the parabolic cylinder
function behind sfm()'s "NG" and "NNAK", which is now computed in R:
a power series in the argument where that does not cancel, and a
peak-centred quadrature otherwise, agreeing with the gsl route to about
1e-10 in log. Below an argument of -sqrt(1400) the gsl code evaluated
its series at a clipped argument, and overflowed at larger shapes; that
behaviour is reproduced exactly rather than corrected, because it is what
keeps the optimizer out of the sigma_v -> 0 corner where the composed
likelihood is unbounded.
psfm() gains the simulated-ML draw controls sfm() has had since 1.2.0
(gap H10): sim_type ("halton", "sobol", "torus", "uniform"),
antithetics, sim_burn, sim_scrambling and sim_prime, for "GTRE"
with estimator = "sml", "TRE", "GTRE_Z" and "TRE_Z". The seed remains
rand.gtre. The panel draws are now built by the same code as the
cross-sectional ones, keeping one block per firm and the shift
randomization. At the defaults the draws and the fits are bitwise identical
to 1.2.0's.
Score-based diagnostics now work for psfm()'s simulated-ML panel models
("GTRE", "TRE", "GTRE_Z", "TRE_Z"). keep_objective = TRUE was
accepted and stored a likelihood, but the closure had no per_obs branch,
so influence_sfa(), vcov(type = "bhhh"), the sandwich methods, TIC()
and vuong() all failed on it. The closure now returns one contribution per
firm, since firms are the independent units, which makes the sandwich
standard errors firm-clustered and the influence diagnostics firm-deletion
diagnostics. bread() scales by the number of firms: scaling by nobs()
with firm-level scores would have overstated the sandwich variance by a
factor of T squared, silently.
data_gen_p() can generate a genuinely non-null first-difference design.
Two new arguments for the y_fd column. fd_draw = "truncated" draws the
firm-level inefficiency from the truncated normal N+(mu, sig_u^2) that Wang
and Ho (2010) and psfm(model_name = "FD") assume; the default, "folded",
keeps the historical |N(mu, sig_u^2)|. The two coincide only at mu = 0.
delta_fd sets the coefficient in u_it = exp(delta_fd * z_it) * u_i*
separately from mu; previously the same mu controlled both, so a design
with delta_z != mu could not be generated. Both defaults reproduce the old
output exactly. Requested by Christopher Parmeter for the convergence study.
?endogeneity_test now carries a size and power table. The test shipped
in 1.2.0 with a theoretical argument for why it should be well calibrated at
the null -- ivsfm() parameterises rho = t/sqrt(1 + t't), so rho = 0 is
an interior point and the delta-method Jacobian is the identity there -- and
a rejection rate measured at a single sample size. Size at one n shows
neither convergence to nominal nor its absence, so the table is a sweep.
Measured over 1,000 replications per cell at
n = 500, 1000, 2000 and 4000: size is 0.043-0.053 against a nominal 5% and
power rises monotonically to one. The variance estimate was checked directly,
not inferred from the rejection rate -- the mean reported standard error of
rho over the actual sampling spread of rho-hat is 0.95 to 1.05 across
every cell. No code changed; the claim is simply now supported.Four approaches to the Waldman (1982) boundary now sit together under one help
topic, ?sfa-wrongskew, which says what each does and when to prefer it. Two
of them are new estimators; the third is a new diagnostic; the fourth
(esfm()) shipped in 1.2.0.
sfm(estimator = "acols") is the alternative corrected OLS of Parmeter
and Zhao (2023, Empirical Economics 64:2831). Corrected OLS with the first
absolute moment E|eps| in place of the third central moment -- the third
moment is the only ingredient that can drive sigma_u to the boundary, so
removing it removes the failure. Implemented for NHN and NE, which are
the two the paper covers.
Replicated against the paper's own Monte Carlo. At n = 200 every COLS
figure matches to 0.005 or better, and the ACOLS Type I failure rates match
too: lambda = 0.2 gives 0.162 against the paper's 0.160, lambda = 0.4
gives 0.151 against 0.146, lambda = 1.5 gives 0.026 against 0.022. The
Normal-Exponential table is reproduced essentially exactly.
One thing the paper leaves open had to be settled. Its recipe is to minimize
the squared moment equation from many starting values, but that equation has
more than one exact root in 18 to 36 per cent of samples, and every root
attains the same objective of zero. sfa keeps the root with the highest
profile likelihood -- a tiebreak, never an optimization -- which reproduces
the paper's failure rates and, unlike "always the first" or "always the
last", is never far off its mean squared errors at either end of the lambda
range.
The NE moment equation is also written through log Phi rather than
2 exp(lam^2/2) Phi(-lam), which overflows above lam = 37. That overflow
is why the paper had to treat sigma_u < 0.03 as zero; here the equation
stays finite to lam = 200.
sfm(estimator = "cmle") is the moment-constrained maximum likelihood of
Zhao and Parmeter (2022, Economics Letters 221:110901). The same
likelihood estimator = "mle" maximizes, subject to two equality
constraints the assumed density implies. Nothing is respecified, and the two
estimators are asymptotically identical; the constraints simply do not hold
at the OLS stationary point, so the search cannot settle there.
The implementation does not need a constrained optimizer. Given beta the
two constraints pin both variance parameters -- in closed form for NHN, by
one scalar root-find for NE -- so the problem is a profile over beta
alone. Clamping at the edge of the admissible set rather than rejecting keeps
that profile continuous, and a fit that ends up clamped is flagged through a
new $boundary component and a warning.
On wrong-skew samples, sigma_u collapses to zero in 95 to 97 per cent of
standard ML fits and 2.5 to 4.5 per cent of constrained ones, against the
paper's claim of a reduction "over 90 per cent". Its beta0 mean-squared
error ratios of 2.07 to 3.40 bracket the paper's 2.15 to 2.76.
skewness_decomp() implements the third-moment decomposition of Bonanno,
De Giovanni and Domma (2017, JPA 47:49), the paper behind copsfm(). It
splits E[(eps - E eps)^3] into the part from the asymmetry of u, the part
from the asymmetry of v, and the part from the dependence between them,
and reports Zenga's median-based measure alongside -- the measure whose sign
says which way the density actually leans.
The point is what it says about the textbook model: a normal v zeroes the
second component and independence zeroes the third, so under the standard
specification the sign of the composed third moment is the sign of
-E[(u - Eu)^3] and nothing else. Run it on an NHN fit and two of the
three components come back identically zero.
The marginals are closed forms; the dependence term is Gauss-Legendre
quadrature on the unit square against the copula density, so it covers every
marginal and copula copsfm() supports rather than only the
exponential/generalized-logistic/FGM combination the paper writes out. It
also settles the two errors 1.2.1 already recorded in those papers' summary
moment formulas: the corrected expressions are now derived and pinned, and
?copsfm gives them.
ttsfm() now reports the three scale parameters on their natural scale,
with names. This is a breaking change to $out and coef(). The rows
are sigma_v, sigma_u and sigma_w. Previously the matrix held raw
optimizer values: the row labelled sigv was log(sigma_v), and in the
homoskedastic case the two one-sided scales were both labelled
(Intercept), so u and w could not be told apart and none of the three
numbers was on the scale its label implied. Code that exponentiated those
rows must stop doing so.
Where sigma_u or sigma_w carries determinants there is no single scale to
report, so those rows stay as coefficients on the link and take Zu. and
Zw. prefixes — which is also what makes the two blocks distinguishable when
they share a covariate name, something the old labelling could not express at
all. Standard errors on the converted rows carry the delta-method factor.
tt_boundary_report() refits a two-tier frontier from several starting
values of sigma_v, scores every solution on the same likelihood, and says
which is highest. Use it when ttsfm() returns sigma_v at essentially zero
and you need to know whether that is the maximum likelihood estimate or an
optimizer that stopped short.
It is a likelihood comparison, not a test, and reports no p-value. The
hypothesis a reader would attach to it sits on the boundary of the parameter
space, where the usual chi-square reference distribution does not apply, and
no appropriate null distribution has been established for this model. The
ordinary linear model on the same data is reported alongside, because a
collapsed sigma_v is not the same claim as "there is no frontier here".
?ttsfm now documents the "TTHN" boundary behaviour: when sigma_u
and sigma_w are close the fit drives sigma_v to zero, and that collapse
is the maximum likelihood estimate. With sigma_u = sigma_w the difference
w - u is already symmetric and heavier-tailed than a normal, so a small
normal v underneath it is very nearly invisible. On a 2,000-observation
draw the profile in sigma_v is flat to 0.000 log units from 1e-6 upward
and decreasing above 0.01, with no interior maximum. Asymmetric
(sigma_u, sigma_w) does not show this. The two-tier structure is
identified; the noise underneath two equal and opposite tiers is not.
sfm(estimator = "cols") returned the frontier coefficients with the wrong
sign on a COST frontier (inefdec = FALSE), in 1.2.0 as released. The
moment path fits on Yc = inefdec_n * Y, which is the production orientation
whatever the caller asked for: a cost frontier y = x'b + v + u is fitted as
-y = x'(-b) + (-v) - u, so what came back was -b. The scale parameters
were right -- v is symmetric and u is the same u -- which is why it went
unnoticed. On a fit with a true intercept of 1 and slope 0.8 it reported
-1.112 and -0.744 where maximum likelihood on the same data gave 1.103 and
0.751. Fixed, and fitted() + residuals() == y is now pinned in both
orientations. Any cost-frontier estimator = "cols" fit made with 1.2.0
should be re-run.npsfm(method = "SZ")'s DEA step now runs from the DEA package rather
than from a linear program built here on lpSolve. DEA moves into
Suggests and lpSolve moves out; the guard and the error message that
names the missing package are unchanged. The two agreed to 6e-12 across
all four returns-to-scale settings and one and two inputs before the swap,
and tests/testthat/test-dea.R now pins the contract sfa needs from DEA
-- orientation, the returns-to-scale spelling, the invariants
method = "SZ" relies on -- rather than re-testing linear algebra that
belongs to DEA.sfm() and psfm() themselves, had titles that were
just their own names, so the manual's index read FinnishElec FinnishElec.
All six now describe what they are.copsfm() gains skewed noise and a choice of inefficiency marginal, via
new vdist and udist arguments (Bonanno and Domma 2022, Mathematics
10:3876). vdist = "glogistic" is a generalized logistic noise carrying its
own skewness parameter alpha_v; vdist = "logistic" fixes it at 1, and
udist = "exponential" replaces the half-normal inefficiency. copula = "independent" is now offered as well, and estimates no dependence parameter
at all, so the skewed-noise frontier can be fitted without a copula.
The point is the wrong-skewness anomaly. The third moment of the composed
error depends on the skew of v and on the dependence between v and u,
not only on the skew of u; where v is symmetric and independent a
positive residual skew has nowhere to go but sigma_u = 0. The paper's four
specifications are four (copula, vdist) pairs -- see ?copsfm.
The paper publishes the composed density in closed form, as four Gauss
hypergeometric terms, for both a production frontier (Theorem 1) and a cost
frontier (Theorem 2). copsfm() computes it by quadrature, and the two agree
to about 1e-14 in log f across alpha_v, theta and both
orientations. Two of the paper's summary moment formulas do not survive the
same check and are not consistent with its own appendix; ?copsfm says which.
sfm(model_name = "NGB2") fits a generalized beta of the second kind
inefficiency (Makiela and Mazur 2022, JPA 58:35-54, section 3), with a
scale and three shape parameters. It nests
most of the rest of the table -- half-normal, exponential, gamma and Weibull
as nu -> Inf limits, and the half-Student t exactly -- so one fit says
which of them the data want. nu is flat above a few hundred, where the
family has reached its generalized-gamma limit, and sfm() warns rather than
reporting the number as an estimate.
The composed density is evaluated by deterministic quadrature with two node sets, one on the inefficiency's scale and one on the noise's, combined by the balance heuristic. It was written on the simulated-ML path first, and that was wrong in a way worth repeating: with 100 draws per proposal the simulated likelihood is biased upward and the optimizer chases the bias, reporting a log-likelihood of -336.9 at parameters whose exact value is -458.7. Four shape parameters give simulation error room to be maximized.
sfm(model_name = "NB", estimator = "cols") fits Carree's (2002)
binomial inefficiency by corrected OLS. It is the only inefficiency
distribution here that can be skewed either way, which is the whole
reason for it: a positive residual skew, which every other model must report
as sigma_u = 0, is here an admissible p > 1/2 and says that most firms
carry considerable inefficiency. The estimator inverts the second, third and
fourth residual moments; it refuses maximum likelihood, and reports Carree's
infeasible region as "no binomial solution" rather than as a number.
copsfm(inefdec = FALSE) fitted the cost frontier with the density
mirrored, in 1.2.0. eps is already sign-normalized, so multiplying the
quadrature node by S as well reflected the composed density: the function
maximized f_{v-u}(y - Xb) where it needed f_{v+u}(y - Xb). On clean cost
data with sigma_u = 1 and sigma_v = 0.4 it returned sigma_u = 0.031 and
sigma_v = 0.69, with the intercept 0.8 too high. The slopes were unaffected,
and so was every production fit (inefdec = TRUE, the default), which is why
the existing tests did not catch it. Any copsfm() cost fit made with
1.2.0 should be re-run.
A parameter on its bound cost copsfm() every standard error in the fit.
The likelihood refused out-of-range draws with .Machine$double.xmax;
optim() differences the objective to build its gradient, and differencing
1.8e308 overflows to a non-finite value, so the final stage aborted and all
standard errors came back NA. It bites whenever a parameter ends up ON a
bound, which a dependence parameter routinely does. Now a large finite
penalty, with the two bounded parameters clamped rather than refused --
matching what sfm()'s NGE/NLN/NW branches already did.
sfm(estimator = "cols") predicted efficiency with the wrong posterior for
"NE" and "NG". E[exp(-u)|eps] is a property of the assumed u, and
the COLS path applied the normal/half-normal posterior regardless of
model_name. Parameter estimates were correct; only exp_u_hat was affected.
It now dispatches, "NG" using the same parabolic-cylinder form its
maximum-likelihood branch uses.
ivsfm(model_name = "IVLIML") could walk into a corner of its own box and
report the result as an ordinary fit. IVLIML is the default estimator, so
this is the one to know about. On the design of Hou, Ramalho and
Roseta-Palma (2025) it returned rho = -0.99503719 with beta2 near 6
against a truth of 0.5 -- on 27 of 300 draws at rho = 0.8 and 6 of 300
at rho = 0.4 -- with no error, no warning and a plausible-looking
sigma_u. "IVCF" failed on 0 of 300 of the same draws.
The cause was the STARTING POINT, not the search. rho = t/sqrt(1 + t't),
and t started at 0 -- which is rho = 0 exactly, the one point where the
frontier and reduced-form blocks are independent and the likelihood carries
no gradient information about which way rho should move. It is also the
centre of the box, because span = pmax(10 * abs(start_v), 10) gives a bound
of +-10 when the start is zero, so the corner t = -10
(rho = -10/sqrt(101)) sat exactly as far from the start as the true
optimum.
IVLIML is now seeded from the IVCF solution. That is not an
approximation borrowed from another model: IVCF is the same likelihood with
Pi and chol(Sigma_xi) pinned at their OLS values, so its answer is a
valid point of IVLIML's own parameter space once those two blocks are
appended. It lands t near its optimum, which both supplies a gradient and
recentres the box away from the corner. Measured over 300 replications at
each of rho = 0, 0.4 and 0.8:
| rho | collapses before | collapses after | RMSE(beta2) before | after | |---|---|---|---|---| | 0.0 | 0/300 | 0/300 | 0.0261 | 0.0268 | | 0.4 | 6/300 | 0/300 | 0.9188 | 0.0263 | | 0.8 | 27/300 | 0/300 | 1.6808 | 0.0241 |
The fit is also about twice as fast, because the optimizer now starts near
the answer. "IVCF" and "C2SLS" are untouched -- the seeding runs only for
IVLIML, and a sub-fit that fails or does not improve on the old start is
discarded rather than used.
The DESCRIPTION now says what the package does. The previous text
("a user-friendly framework ... extensive flexibility in specification and
estimation techniques") named no model, no method and no reference. It now
names the model families and cites Aigner, Lovell and Schmidt (1977),
Meeusen and van den Broeck (1977), Greene (2005), Colombi, Kumbhakar,
Martini and Vittadini (2014) and Amsler, Prokhorov and Schmidt (2016) with
DOIs.
A feature release. In brief:
A Wald test for whether the endogeneity correction was needed,
endogeneity_test(), from Hou, Ramalho and Roseta-Palma (2025). On an
ivsfm() fit it tests H0: rho = 0 -- that the noise is uncorrelated with
the reduced-form errors, so a plain sfm() fit would have been consistent.
The null is INTERIOR, since rho lives in the open unit ball, so this is an
ordinary chi-square and not the chi-bar-square mixture
inefficiency_test() needs. Size 0.046 at a nominal 5% over 1000
replications; power 1.000 at rho = 0.4 and 0.8.
ivsfm() now reports standard errors for rho. They were NA, on the
ground that rho = t/sqrt(1 + t't) has a non-diagonal Jacobian and a
delta-method value computed as though it were diagonal would be wrong. The
first half was right and the conclusion was not: the Jacobian is
J = (I - rho rho')/s, so Var(rho) = J Var(t) J' exactly. Checked against
a numerical Jacobian, and validated end-to-end by the size of the test above
-- a wrong Jacobian would show up there as a mis-sized test. The full
covariance travels with the fit as $vcov_rho, which the joint test needs.
ttsfm(model_name = "TTHN") no longer fabricates a log-density it cannot
compute. The likelihood needs D = Phi2(x1, 0; rho1) - Phi2(x2, 0; rho2),
a difference of two bivariate normal CDFs of the same order of magnitude, and
in parts of the parameter space it cancels completely: at
sigma_v = 0.3, sigma_u = 1, sigma_w = 0.2 on a 400-observation draw one
observation's D comes back as exactly 0 while both CDFs are O(0.1).
pmax(D, .Machine$double.xmin) did not repair that -- it invented
log(D) = -708.4 for the observation, and that single observation moved the
summed objective by 715 log-units in a surface whose real curvature is a
few units per 0.1 step in log sigma. A cliff that size is what L-BFGS-B
reports as ABNORMAL_TERMINATION_IN_LNSRCH. D is now rejected when it
falls below the rounding error of its own subtraction, which leaves genuine
small probabilities alone: a D of 6e-94 where both CDFs are also ~1e-94 has
full significance and is kept.
Three smaller repairs in the same branch: the exponentials for
sigma_v/sigma_u/sigma_w are clipped at EXP_CLIP_UPPER as the TTNE
branch already clipped its own; rho is held strictly inside (-1, 1),
which sigma_v small enough does reach; and the substitution of
-sqrt(.Machine$double.xmax / n) -- about -9.5e152 at n = 200, 140 orders
of magnitude past the finite penalty it sat beside -- is gone, replaced by
that penalty. It had also mapped a +Inf log-density to a large NEGATIVE
one, so an unusable draw could be rewarded.
Measured on eight matched fits: one that previously failed outright now
converges, and one that previously returned sigma_w = 0.044 against a truth
of 1 now returns 1.047 with a better likelihood. The other six are
unchanged. sigma_v still collapses toward zero on some draws; that is a
separate problem and is not fixed.
ttsfm() no longer throws away a fit because of the optimizer's exit
code. It was the only entry point calling stop() on ANY non-zero
convergence code, including L-BFGS-B's code 52,
ABNORMAL_TERMINATION_IN_LNSRCH, which it returns whenever its line search
meets a discontinuity and usually having improved the objective. The guard
was written when nothing checked opt$convergence at all; opt.optim() now
rebuilds at the stage-2 point whenever the value or the Hessian is
non-finite, so what was left of it only destroyed usable fits. Two of eight
TTHN test fits returned code 52 with sound estimates where the stop() had
aborted them.
What skewness and kurtosis a pair of distributions can actually produce,
moment_range(), from Papadopoulos and Parmeter (2021). spec_test() already
shipped their formal statistic; this is the range check that comes before it.
Because the variance share R lies in (0,1), the composed error can never
be more skewed than u itself -- Normal-Half-Normal is confined to
gamma1 in (-0.995, 0), Normal-Exponential to (-2, 0) -- so residual
moments outside that set refute the pair with no test statistic involved.
That is their Type II failure, alongside the Type I failure of wrong skew.
Five noise and five inefficiency distributions, including the two whose
shape is a free parameter: the truncated normal sweeps
gamma1(u) in (0, 2) and gamma2(u) in (-0.243, 6) and is the only u
here whose composed error can be platykurtic under normal noise, while the
gamma is unbounded in both and so can never suffer a Type II failure.
The kurtosis range is not the interval between the two components: with
leptokurtic noise it dips below both, to 0.674 for Laplace noise with a
half-normal u against the half-normal's own 0.869.
It is a diagnostic, not a test, and the print method says so. PP2021's
Tables 1-5 replicate -- their population rows exactly, and Tables 2, 3 and 5
cell for cell -- and they show the range check over-rejects: with correctly
specified data at n = 200 and SNR 3, the sample excess kurtosis leaves its
theoretical range 26% of the time. Use spec_test() for the p-value.
A frontier that tolerates the wrong skewness, esfm(), from Hafner,
Manner and Simar (2018). When the sample skewness of the residuals comes out
positive -- a small-sample accident, not evidence against the model -- the
classical MLE collapses: sigma_u is exactly zero, every firm is fully
efficient and the fit says nothing. This model keeps the frontier and widens
the distribution of u, with one parameter gamma whose SIZE is the scale
of inefficiency and whose SIGN chooses the direction of skewness. The
classical model is nested at gamma > 0.
On a 50-observation wrongly skewed draw, sfm(model_name = "NHN") returns
sigma_u = 0.0003 and mean efficiency 0.99976; esfm() returns
gamma = -0.57 with the slopes still near truth.
Because gamma = 0 is an interior point rather than a boundary, the LR
test of "no inefficiency" is an ordinary chi2(1) -- not the chi-bar-square
mixture inefficiency_test() needs for the classical model. Size over 2000
replications on the paper's design, at a nominal 5%: 0.178 / 0.082 / 0.054 /
0.061 / 0.058 for n = 50 / 100 / 200 / 500 / 1000, against the paper's
0.152 / 0.076 / 0.061 / 0.054 / 0.053. Badly oversized below n = 100, usable
from n = 200 -- their own conclusion.
It is not a free improvement everywhere. The paper reports smaller bias
in mean efficiency than the classical model even when the population
skewness has the correct sign, and that reproduces when inefficiency is not
small: at gamma = 0.5 the bias is -0.0001 at n = 100 against the classical
+0.021. But at gamma = 0.3, where 22% of samples are wrongly skewed, this
model overshoots the other way (-0.040 against +0.036) and is the worse of
the two at n = 50 (-0.080 against +0.044). A small true gamma lets it fit
sizeable inefficiency of either sign to noise. Use it when the wrong
skewness is the problem.
Kim and Schmidt (2008) two-step test, uhet_test(). The procedure
applied work actually runs -- fit the frontier, regress the JLMS predictor
on firm characteristics, test the slopes -- is invalid as usually performed,
because u-hat is a generated dependent variable and the first-step
estimation error enters the second-step variance. The correction is not
distribution-free, so "NHN" and "NE" are implemented, the two cases the
paper works through.
What it buys: with z independent of x the corrected and naive tests
agree, as the theory requires. With z strongly correlated with x
(rho = 0.9) the naive standard error is 52% too large at n = 800 and
the naive test's size collapses to 0.000 at a nominal 5%, while the
corrected standard error matches the true sampling spread to within half a
percent (ratio 0.995) and keeps size 0.032.
Caveat, stated plainly: the corrected test is conservative in these
simulations. On the paper's own base case it rejects 0.024/0.036/0.028 at
N = 200/500/1000 against a nominal 0.05, where the paper reports 0.0495 at
N = 200. The naive branch does reproduce the paper's figure for the
uncorrected test (0.023-0.029 against their 0.0237). The implementation was
checked against equations (5)-(13), the first-step scores agree with an
analytic half-normal score to 1e-11, gamma-hat is unbiased and Gaussian,
and an oracle test using the true sampling standard deviation has correct
size (0.045) -- so the gap is in the variance estimate, not in gamma-hat
or the reference distribution. Two candidate explanations were tested and
both refuted: the homoskedastic second-step variance (the paper's "OLS"
against its "OLS-H") changes nothing, and neither does using the Hessian
information rather than the outer product when forming r_i. The gap is
unexplained rather than merely unverified. Read a rejection as meaningful
and a non-rejection as weak evidence.
Fits on the wrong-skew boundary are refused rather than answered: every
E[u|eps] is then identical, so there is no estimated inefficiency for z
to explain, and the information matrix is singular there (Waldman 1982).
Chen and Wang (2012) centered-residuals moment test, cw_test(). Unlike
gof_test(), it never evaluates the composed-error density or distribution
function: it compares the empirical characteristic function of the centred
residuals with its theoretical value, which follows from the characteristic
functions of the two components separately -- in closed form for many pairs
whose convolution is not, the normal-gamma model being the standard example.
Centring is what makes it work. A frontier's intercept is not identified
separately from E[u], so the least-squares intercept estimates
alpha - E[u]; centring cancels it, which is why the test is valid off an
ordinary regression rather than requiring maximum likelihood. The variance
corrects for both nuisances -- estimating the scales, and the centring
itself. Omitting the second understates it and oversizes the test.
Size at a nominal 5% over 500 replications, half-normal data:
| test | n | tau=1 | tau=1.5 | tau=2 | |---|---|---|---|---| | cosine | 500 | 0.070 | 0.054 | 0.050 | | cosine | 2000 | 0.046 | 0.044 | 0.050 | | sine | 500 | 0.092 | 0.064 | 0.042 | | sine | 2000 | 0.068 | 0.062 | 0.064 |
Power against exponential inefficiency at tau = 1 is 0.231 (n = 500) and
0.952 (n = 2000) for the cosine test against 0.108 and 0.706 for the sine
test. Both of the authors' conclusions reproduce: the cosine test is the
less sensitive to tau and the more powerful, so it is the default.
tau defaults to a single frequency: passing several at once gives
rejection rates of 0.09-0.15 at a nominal 5%, because the cosine moments at
nearby frequencies are near-collinear and the covariance matrix is then
near-singular. The authors report the same for their combined tau.
A wrongly signed third moment leaves the moment equations with no admissible solution, and that is an error rather than a number.
The Gauss-Legendre rule is now cached. It depends on nothing but the
node count and was rebuilt by an eigendecomposition on every call --
hundreds of times per cw_test(), since numDeriv differences a
quadrature-based moment function. Caching also keeps eigen() out of forked
workers: macOS's threaded BLAS is not fork-safe and mclapply() segfaulted
inside it before this was added.
Goodness-of-fit tests for the assumed inefficiency distribution,
gof_test(), from Wang, Amsler and Schmidt (2011). The distribution of u
is the assumption in this model least often defended and least often tested,
and it is testable: hold the normality of v fixed and the assumed u
implies a distribution for the composed error, so rejecting that is
rejecting the assumed u. Kolmogorov-Smirnov and Pearson chi-square, both
on the composed error rather than on u-hat -- which is the same test,
since u-hat = E[u|eps] is monotonic in eps, but far easier, and it
avoids the trap the paper warns about: comparing the spread of u-hat
with the assumed density of u is a mistake, not a diagnostic, because
those are different distributions.
Parameter estimation is not ignored. The default parametric bootstrap copies
the estimation step exactly -- each replication draws from the fitted model,
rebuilds the response, refits, forms its own residuals and recomputes the
statistic. The "asymptotic" option gives a chi-square p-value that the
authors note is conservative at the MLE, and returns NA rather than a
KS p-value, because with parameters estimated the Kolmogorov distribution
does not apply and Bai's (2003) transformation is not implemented.
Size and power over 200 replications (half-normal data for size, exponential data fitted as half-normal for power), at nominal 0.10 / 0.05:
| | n | KS | chi2 | |---|---|---|---| | size | 200 | 0.065 / 0.025 | 0.085 / 0.065 | | size | 800 | 0.070 / 0.040 | 0.065 / 0.010 | | power | 200 | 0.505 / 0.395 | 0.235 / 0.145 | | power | 800 | 0.985 / 0.965 | 0.810 / 0.655 | | power | 2000 | 1.000 / 1.000 | 0.990 / 0.985 |
Both hold their size, mildly conservatively, and KS dominates chi-square at every sample size -- the paper's own conclusion, and the reason it is the better default.
Coelli (1995)'s hypothesis tests, and his COLS standard errors. New
inefficiency_test() tests H0: no technical inefficiency, reporting the
one-sided LR test, the naive LR test, the Wald ratio and the third-moment
test side by side. The point of the paper is that the two tests most often
reported are the two with the wrong size: H0 puts gamma on the boundary, so
the LR statistic follows the Gourieroux-Holly-Monfort mixture
0.5 chi2(0) + 0.5 chi2(1) (5% critical value 2.71, not 3.84) and the Wald
ratio is not asymptotically normal at all. Measured here over 400
replications per cell at a nominal 5%, with sigma_u = 0:
| n | LR one-sided | LR naive | Wald | M3T | |---|---|---|---|---| | 100 | 0.060 | 0.030 | 0.230 | 0.055 | | 400 | 0.055 | 0.033 | 0.217 | 0.045 | | 800 | 0.060 | 0.035 | 0.182 | 0.062 |
The one-sided LR and third-moment tests hold their size; the naive LR test
rejects about half as often as it should; the Wald test rejects three to four
times too often and does not improve with n. Under the alternative the
one-sided LR test has the better power of the two valid tests at every sample
size, which is Coelli's recommendation. NTN fits get the two-restriction
mixture 0.25/0.5/0.25 over chi2(0,1,2).
sfm(estimator = "cols") now reports Coelli's analytic standard errors.
It previously returned NA for sigv and sigu unless a bootstrap was
requested, because the OLS errors are not appropriate for them: those
parameters are non-linear functions of the residual moments m2 and m3,
not regression coefficients. Coelli's Appendix 1 (A15-A16) gives the correct
delta-method variance and it is now filled in for NHN; the intercept
additionally carries the error in its E[u] shift. NE and NG keep NA
and the bootstrap, because the moment inversion -- and therefore the
variance -- is distribution-specific.
Checked against the nonparametric bootstrap that was already there, on the
same data: the analytic/bootstrap ratio for sigu is 1.047 at n = 200,
1.011 at n = 800 and 0.984 at n = 3200, which is an asymptotic
approximation behaving as one should.
Coelli's A13 inversion turned out to be the moment inversion the package
already carried as .gtre_two_step() (equal to ten decimal places, not
merely similar), and his A15 variance is the delta method whose GTRE_SEQ2
implementation is corrected elsewhere in this release -- including the
factor of 2 on the covariance term, which A15 independently confirms, and
A16's mu6 - mu3^2 - 6 mu2 mu4 + 9 mu2^3, which confirms the other two
repairs. The central moments now live in one shared helper,
.nhn_central_moments(), rather than being derived twice.
A composed-error CDF for every cross-sectional model, pcomposed_model()
and composed_cdf(). pcomposed() covers the half-normal only, which is
enough for copulas but not for goodness-of-fit testing, where the CDF is
needed for whichever inefficiency distribution the model assumes. All
thirteen are covered: NHN, NTN, NE, NR, NU, NGE, NLN, NW, NG, NNAK, TSL,
THT and tHN.
Three parameterizations do not read the way their names suggest and are
taken from each likelihood rather than from the label: THT lists sigu
before sigv, the only model that does; NG's sigu is the gamma
scale and its mu the shape; NNAK's sigu is the Nakagami spread,
so Omega = sigu^2; and NLN's mu is a meanlog.
THT is not an independent convolution. Tancredi's composed error is
skew-t, a scale mixture in which v and u are divided by the same
sqrt(V/a). Treating it as t noise plus an independent half-normal -- which
is exactly what tHN is -- gets the log density wrong by up to 1.74. The
first draft did precisely that, and the check below is what caught it.
Verified three ways: against pcomposed() for the half-normal, agreeing to
3.5e-15 relative including a lower tail at log F = -209; against a
four-million-draw simulation for all thirteen models, maximum absolute error
4e-4 against a Monte Carlo standard error of 7.5e-4; and by numerically
differentiating the CDF and comparing with each model's own stored
log-density, which agrees to about 1e-8 for eleven of them. NLN and NW
differ by ~3e-3 because their likelihoods are simulated, not because the
models disagree.
psfm() died with "system is exactly singular" when a factor was only
partly collinear between individuals. Reported from use: factor(year) on
an unbalanced panel, where the between-individual means of the year dummies
are rank deficient. The guard added in 1.1.5 detected this correctly and
printed an accurate warning, then did nothing, because it could only express
the fix by deleting whole formula terms -- and 11 of factor(year)'s 24
dummies being collinear names no term that can be deleted without discarding
the 13 identified ones too. plm::ercomp() was then handed the same singular
design and failed inside solve(). The starting-value design is now reduced
at column granularity (.re_start_design()), so the identified dummies
survive, the dropped ones take their starting values from a pooled OLS fit as
documented, and the likelihood still uses the formula as written. This was a
no-op for every partly-collinear factor, which is the common case -- a fully
collinear term was the only shape the old path could handle.
Two consequences of that path are fixed with it. collinear_action = "warn_drop" silently dropped nothing in the same situation, and now says so
instead. And "GTRE_SEQ1"/"GTRE_SEQ2" took the frontier standard errors
from the starting-value fit by position, which errored on the length
mismatch once the design was reduced; they are matched by name, with NA for
columns that fit could not identify rather than a neighbour's value.
"GTRE_SEQ1" and "GTRE_SEQ2" are inconsistent for fixed T, and now say
so. Both hand plm's random-effects output to a second stage that treats it
as draws from the latent composite errors. It is not: alpha_hat is a
shrunken BLUP and eps_hat a quasi-demeaned residual, which attenuates the
third central moment by about 0.81 at \eqn{T = 10} on the transient side and
scales the persistent side by \eqn{c^2} and \eqn{c^3}. No factor tends to 1 as
\eqn{N} grows at fixed \eqn{T}, so both converge to the wrong constants --
\eqn{\sigma^2_{hr}} settles near 0.68 against a true 0.80. Feeding the same two
second stages the latent draws instead recovers the truth at the \eqn{\sqrt{N}}
rate, which puts the fault in the shared first stage rather than in the
likelihood or the moment inversion. Each call now warns and points at
model_name = "GTRE". This is a documentation and diagnosis change: the
estimators are unchanged, and are the large-T procedures they always were.
"GTRE_SEQ2"'s standard errors were wrong in four ways. The variance of
the third central moment used mu_3^3 for \eqn{\mu_3^2} and 9 * mu_2 for
\eqn{9\mu_2^3} (also dimensionally inconsistent); the delta-method
cross-terms were missing their factor of 2; beta_0's error added a standard
error to a sum of variances instead of its square; and both
\eqn{\partial\gamma} derivatives were garbled, squaring neither denominator and
carrying \eqn{m_3} where \eqn{m_2} belongs. Against the Monte Carlo sampling
standard deviation the reported error for gamma_uv was 16 times too
large and sigmaSq_uv 1.37 times; corrected, both land within 2%. The
algebra now lives in .gtre_two_step_se() so it is unit-tested, and returns
NA at the wrong-skew boundary where it is not differentiable rather than
Inf. "GTRE_SEQ2" also reported a lambda computed from a different
(likelihood) decomposition than the sigmas printed beside it; it now comes
from its own.
ttsfm() returned an infinite penalty on unusable parameter draws, in
seven places. TTNE, TTHN and TTNLS each returned
.Machine$double.xmax when a draw produced non-finite contributions. That is
the construct sfm() was converted away from earlier in this release after
"NGE" died on 3 of 45 fits: optim() differences the objective to get its
gradient, and differencing 1.8e308 overflows to a non-finite value, which does
not steer the search away -- it aborts the fit. TTHN was failing at
\eqn{n = 400} with ABNORMAL_TERMINATION_IN_LNSRCH, which is what L-BFGS-B
does when handed gradient information it cannot use. All seven now return a
finite penalty.
The regression test that exists to catch exactly this construct only ever read
sfm.R, which is how ttsfm.R kept all seven for as long as the test
existed. It now reads every file carrying a likelihood closure.
ttsfm() no longer calls the bivariate normal CDF once per observation.
The TTHN density was evaluated with mapply() over observations, so one
mnormt::pmnorm() call per observation, twice per likelihood evaluation --
and with a numerically differenced gradient, \eqn{2n(2p+1)} calls per
optimiser iteration. The correlation argument varies only through
\eqn{\sigma_u} and \eqn{\sigma_w}, so under homoskedasticity it is one value
repeated and the whole vector is a single call. Grouped accordingly, with a
fallback to the per-observation form when the correlation genuinely varies
per observation -- there, grouping is slower, measured at 0.60x. Results are
bit-identical to the previous form. The gain is about 1.4x, not more:
pmnorm() still loops internally, so the call count falls but the work does
not. TTNLS separately computed each of its two bivariate terms twice; those
are now computed once.
selsfm() gains a second model, model_name = "kts": Kumbhakar, Tsionas
and Sipilainen (2009), joint estimation of technology choice and technical
efficiency. Not a variant of Greene (2010), which this entry point already
fitted, but a different model. Greene has one frontier observed only for the
selected, with selection correlated with the noise. Here both technologies
are observed, each with its own frontier and its own two scales, and the
choice depends on inefficiency itself -- u lowers output through the
frontier and shifts the adoption decision through the choice equation.
Neither two-step order works, which is the paper's own argument for
single-step maximum likelihood: the choice equation cannot be a probit
because u is unobserved, and fitting the frontiers first ignores the
endogeneity of the choice. Because u's distribution depends on the regime
while the regime's probability depends on u, the marginal density of u is
a mixture of the two half-normals and the observation's density is the same
mixture of two one-dimensional integrals. n_nodes controls the quadrature.
copsfm() goes from two copula families to fifteen: Frank, Clayton,
Gumbel and Joe join Gaussian and FGM, each of the last three with 90/180/270
rotations, since Clayton, Gumbel and Joe carry only positive dependence and
nothing rules out a negative association between noise and inefficiency.
Every density is verified against the second mixed partial of its own CDF, it
integrates to 1 over the unit square, and it is exactly 1 at the independence
parameter.
copsfm() now warns when you choose a family that does not recover its own
dependence parameter, and most of them do not. Fitting each family to 25
samples generated from itself at n = 400, the share of fits returning the
independence boundary was 0% for Frank and Clayton but 36% for Gumbel, 40%
for Joe, 56% for Clayton rotated 270 and 60% for Gumbel rotated 90. On data
from a Gumbel copula with Spearman 0.685 at n = 2000, every family
including the true one returns the boundary and their log-likelihoods differ
by less than 0.04: the likelihood is close to flat in the dependence
parameter. That is a property of the model, not of the implementation. The
warning quotes the measured rate; the rotations that were not measured say
so. Prefer "frank" or "clayton".
zsfm(model_name = "ZISF") could be trapped at the boundary where the mass
point vanishes, returning a much worse fit without saying so. Since
P(fully efficient) = exp(-gamma), gamma -> Inf is the boundary at which
the model collapses to an ordinary frontier, and it carries its own local
optimum. On the package's convergence design the default start fell into it
on roughly 1 replication in 5000: one fit returned gamma = 9.93
(p = 5e-5) with a log-likelihood 188.7 points below what a start from
anywhere in the interior reaches. Nothing in the returned object flagged it.
zsfm() now tries several starting values for gamma and keeps the best
attained objective, so it can only improve a fit. The guard applies to
"ZISF" only: "ZISF_Z" has no scalar gamma -- its layout is
(sigv, sigu, beta, z-block) with the mixing parameterised through the z
coefficients -- so the same perturbation there would move the noise scale.
Worth knowing if you have compared this model's convergence behaviour before:
that single fit inflated its sample size's mean squared error 144-fold and
flattened the log-MSE-on-log-n slope for gamma to zero, while four of the
five sample sizes were textbook (MSE 0.0047, 0.0019, 0.00095, --, 0.00065).
The sigma_u -> 0 repair listed further down opened a second catastrophic
cancellation, at sigma_v -> 0, which is also fixed here. To be clear about
the scope, since it affects whether anyone need re-check old results: this
defect was introduced by that repair and existed only in this release's own
development. No published version of the package has it -- 1.1.5 and
earlier have the sigma_u defect instead, and .log_phi_tilt() does not
exist in them. It is recorded because the mechanism is worth knowing and
because the fix changes NE's numerics.
That repair works through .log_phi_tilt(), which returns log Phi(z) + z^2/2. For an observation
with eps < 0 and sigma_v small, z goes large positive, log Phi(z)
is ~0, and so the returned value is essentially z^2/2 -- which the NE
branch then subtracted eps^2/(2 sigma_v^2) from. At sigma_v = 1e-8 both
are about 5e15, where consecutive doubles are ~1 apart, so the result was
rounding noise: the summed negative log-likelihood read -5188 where the
true value is +4398, an apparently excellent fit in a narrow band that a
minimiser walks straight into. The domain guard did not catch it, firing only
at sigma_v <= 0 while the band sits at 1e-9 to 1e-7.
The consequence was not a visible failure. nlminb found the band, its
"only accept an improvement" guard passed the point through because the
objective had improved, and bobyqa and optim could not escape the corner;
optim's bounds, rebuilt around the corrupted point, then pinned the
frontier coefficients. Such fits came back as ordinary sfareg objects with
convergence = 0 -- correctly, since a bound-constrained stationary point is
a legitimate success -- while carrying a gradient in the hundreds.
The subtraction is analytic: z^2/2 = eps^2/(2 sv^2) + eps/su + sv^2/(2 su^2)
identically, so above NE_TILT_SWITCH the expression collapses to the 1.1.5
tilt with no large intermediate formed. The two forms agree to 1e-13 across
the switch, and below it the tilt form -- and hence the sigma_u -> 0 fix --
is untouched. Over 500 replications at lambda = 3.33, n = 1000, the rate of
fits landing more than 0.3 from the truth falls from 3.8% to 1.0%, the
log-likelihood improves on 18 and worsens on none.
The root-n convergence sweep was then re-run in full, 1000 replications at
each of n = 1000 to 5000, and "NE" now passes on the mean basis, which is
what it used to fail: slopes -1.023 (sigma_v), -1.044 (sigma_u), -1.025,
-1.033 and -1.009 (the frontier coefficients), every R-squared at least 0.989,
and no failures in 5000 fits. Before the fix those slopes were -0.488, +0.262
-- positive, i.e. error growing with the sample -- -0.745, -0.443 and -0.165.
This also resolves a dilemma the convergence notes had recorded and left open.
They found that reverting the likelihood to its 1.1.5 form fixed the tail but
would reintroduce the sigma_u -> 0 defect, and concluded "do not revert, fix
the interaction instead". Above the switch this fix is the 1.1.5 algebra, so
the re-run reproduces that configuration's slopes digit for digit while
keeping the sigma_u protection. Both boundaries are now covered.
"NGE" shares .log_phi_tilt() and was checked for the same defect; it does
not have it.
Robust estimation gains its missing half. sfm(robust = ) already fitted
the Normal--Half-Normal frontier by maximum \eqn{L_q}-likelihood, the power-Psi
criterion and minimum density-power divergence, but the tuning parameter had
to be supplied by hand with no guidance and no diagnostic. Four new functions
close that, from Bernstein, Parmeter and Wright (2026):
hscore_select() chooses the tuning parameter by minimising the Hyvarinen
score (Sugasawa and Yonekura, 2021) over a coarse-to-fine grid that includes
the maximum-likelihood endpoint, so the criterion is free to decline
robustification. It guards two silent failure modes: candidates lost to
numerical overflow, and warm-start capture in the degenerate
sigma_u -> 0 basin.hscore() evaluates the criterion, in logarithms. The natural-scale
expression raises the fitted density to the power c - 2, so one
observation whose density underflows makes the score non-finite; because the
surviving candidates are those closest to maximum likelihood, the failure
biases selection toward too little robustness without any warning.
stable = FALSE reproduces the naive evaluation for comparison.calibrate_c() gives the fixed alternative: the tuning value at which a
residual three scale units out retains a stated share of its clean-data
influence. The influence ratio is not monotone for the
Fisher-consistency-corrected criteria -- Psi and MDPD cross the target twice
-- so every root is found and reported in a "roots" attribute rather than
resolved silently.density_weights() returns the per-observation weight the estimator
applies, documented with what it can and cannot detect: it responds to
response contamination and not to a mis-recorded regressor, which can be
highly influential and still receive a weight near one.The power-Psi and MDPD objectives satisfy obj_MDPD = (1 + c) * obj_Psi
identically, so they share a maximiser and are the same estimator at a common
c. Both remain available as peer methods and calibrate_c() returns the
same value for either.
influence_sfa(), an influence-function diagnostic for any "sfareg"
fit made with keep_objective = TRUE, after Stead, Wheat and Greene (2023).
The standard outlier rules do not transfer to this model -- they ignore the
asymmetry of the composed error, so a genuinely inefficient firm reads as an
outlier and is discarded. It reports the empirical influence function, a
Cook's-distance-shaped case influence, and two sensitivities. Compare
models on sensitivity_std, the self-standardised one: the raw sup-norm
depends on the parameterisation, and on a clean sample reads 72.6 for "NHN"
against 1322.5 for "tHN", which says nothing about robustness. On the
standardised scale, contaminating one response takes "NHN" from 13.5 to 32.3
and "tHN" from 13.1 to 17.2.
A fifth model-fitting entry point, lcsfm() — the latent class
stochastic frontier, with n_class coexisting technologies.
New accessors. efficiency() returns technical efficiency with a choice
of predictor ("bc", "jlms", "mode") and of scale for the dependent
variable (logDepVar); the mode is the only one that reaches exactly 1, which
it does wherever the posterior mean of inefficiency is negative.
meanefficiency() returns the model-implied E[exp(-U)] for ten
distributions -- seven in closed form -- together with mean efficiency among
the most efficient p, which is a property of the fitted model rather than of
the sample. simulation_se() reports how much of a simulated-ML standard
error is simulation noise rather than sampling noise.
sfm() gains four arguments. weights/wscale for observation weights
(a weight of 2 gives the same estimates as including the row twice);
start_from to seed a hard model from a simpler fitted one, matched by
parameter name; scaling for the Wang-Schmidt scaling-property model on
"NTN"; and shapehet, an experimental covariate-dependent shape for
"NG"/"NNAK" whose limitations are set out in ?sfm.
vcov(fit, type = "bhhh") returns the outer-product-of-gradients
covariance, which is defined even where the Hessian is not.
Regression tables. texreg::screenreg(), texreg() and htmlreg() now
render "sfareg" fits directly.
Two fixes to things that were shipped but did not work. The sandwich
methods were defined without being registered, so sandwich::estfun(),
vcovCL() and lmtest::coeftest() all failed on an installed package --
every test passed because devtools::load_all() exports everything. And
sim_seed silently did nothing under the default sim_type = "halton",
leaving the cross-sectional simulated-ML models with unrandomized draws;
they are now shifted like the panel models. Fixing the second is what makes
simulation_se() meaningful.
"NNAK" no longer fails outright at larger shape values.
gsl::hyperg_U throws above order ~16 and the vectorized call meant one bad
element nulled the whole vector. There is now a fallback that evaluates the
parabolic cylinder function by its integral representation, agreeing with the
route it replaces to 1e-12 where both work.
A seventh entry point, ivsfm() — stochastic frontier estimation when a
regressor is endogenous, in the sense of being correlated with the noise.
Three estimators from Amsler, Prokhorov and Schmidt (2016): "IVLIML" (full
maximum likelihood over the frontier and the reduced form jointly), "IVCF"
(the two-step control function of Kutlu 2010) and "C2SLS" (corrected 2SLS,
which is moment-based and so reports no log-likelihood). Supplying uhet
gives the model of Amsler, Prokhorov and Schmidt (2017), in which
environmental variables entering the inefficiency scale may themselves be
endogenous. Ignoring endogeneity is not a small matter: on simulated data
with a correlation of 0.6 between the noise and a regressor's reduced-form
error, an ordinary sfm() fit carries a bias of 0.133 in that regressor's
coefficient where ivsfm() carries 0.002. Note that "IVCF"'s conventional
standard errors are understated unless the correlation is zero, and that the
case where the regressor is correlated with the inefficiency rather than
the noise is not covered.
A sixth entry point, selsfm() — Greene's (2010) stochastic frontier
with a correction for sample selection, for the case where the units in the
sample are there because of something correlated with the frontier's own
noise. The selection and frontier equations are supplied as separate
selection and frontier arguments. Estimation is Greene's two-step
scheme: probit on the full sample, then maximum simulated likelihood on the
selected subsample, with efficiency scores from the simulated Bayes rule
over the same draws. Note that the standard errors are conditional on the
first-stage estimate and do not carry the Murphy–Topel correction; bootstrap
over both stages when inference on rho matters.
Six new panel estimators, all classical rather than maximum likelihood
and none assuming a distribution for inefficiency: "CSS"
(Cornwell-Schmidt-Sickles 1990), "LS" (Lee-Schmidt 1993), "KSS"
(Kneip-Sickles-Song 2012), and "SSRE"/"SSCRE", which complete the
Schmidt-Sickles family that "SSFE" already began.
Heteroskedasticity in more than one error component, as named formula
arguments to sfm(): vhet, uhet and muhet. The last, with
model_name = "NTN", is Battese-Coelli (1995).
Three corrections that change numeric output: "NE"/"NGE" could
return a positive log-likelihood with sigma_u at zero; nobs() counted
rows supplied rather than rows used, so BIC() used the wrong n; and
"NGE" aborted outright on about 7% of small samples.
Plus z_link for sfm() and ttsfm(), marginal_effects() for panel
fits, efficiency_ci(), and sfm(model_name = "TSL").
Nothing in this release changes the meaning of an existing argument or
model_name. Every addition is a new argument defaulting to the previous
behaviour, or a new model_name value.
Heteroskedasticity in more than one component: new vhet, uhet and
muhet arguments to sfm(). sigma_v, sigma_u and (for "NTN") the
pre-truncation mean of u may now each be driven by their own covariates,
in any combination:
sfm(y ~ x1 + x2 | z_u, vhet = ~ z_v, model_name = "NHN_Z")
sfm(y ~ x1 + x2, muhet = ~ z_mu, model_name = "NTN") # Battese-Coelli 1995
Named formulas rather than further pipe segments, because pipe POSITION
already means different things in different families -- the second segment
is sigma_w in ttsfm() but sigma_h in psfm()'s GTRE_Z -- and a
fourth position would not be readable. The existing | z keeps its meaning,
so nothing that worked before changes. uhet is the same specification as
| z and refuses if both are given; it exists so that "NTN", which takes
no pipe segment, can still have a heteroskedastic sigma_u.
Available for NHN, NHN_Z, NE, NE_Z and NTN. muhet requires
NTN, the only family in which a pre-truncation mean exists. Recovery on
simulated data with all three blocks non-constant, n = 1500: every one of
the six parameters within one standard error of truth, and the two z-links
agree to 1e-6 in log-likelihood with deltas differing by exactly the factor
of two the reparameterization implies.
Two reasons this matters beyond completeness. Ignored heteroskedasticity in
v biases the estimated FRONTIER, not just its precision (Caudill, Ford and
Gropper 1995), because the composed error's mean depends on the scales. And
vhet is one of the few specifications that can absorb apparent wrong
skewness without forcing sigma_u onto the zero boundary.
These fits carry every positive quantity on the log scale, so the parameter
vector is unconstrained and no variance can be driven onto a boundary by the
optimizer. marginal_effects() reads them with no extra arguments.
estimator = "cols" and the robust divergences are moment-based and
homoskedastic respectively; both refuse rather than silently ignoring the
new arguments.
print() and summary() now say when a variance scale is on the zero
boundary, for sfm()'s sigma_u under wrong skewness, for "tHN"'s
sigma_u absorbed by heavy-tailed noise, and for psfm()'s two persistent
scales. Each says which of those it is, since they mean different things.
The condition was already warned about at fit time, but a warning is a
one-off: it does not survive saveRDS(), a fresh session, or a loop that
suppressed warnings to keep the console readable. What people look at
afterwards is the printed table, and there a collapsed scale appears as an
ordinary coefficient of 0.0009 with nothing to say that it is a boundary
solution, that the split it belongs to is unidentified in this sample, or
that the intercept moved to absorb the difference.
psfm_bootstrap() warns when asked to bootstrap a boundary fit. A
parametric bootstrap resamples from the fitted model, so a fit whose
persistent scale has collapsed resamples from a data-generating process in
which that component is absent -- and the bootstrap is inconsistent on the
boundary of the parameter space in any case. The intervals would understate
the uncertainty in the persistent split rather than represent it, while
looking perfectly ordinary.
Since one of GTRE's two persistent scales collapses in roughly a third of samples, this is expected to fire in ordinary use. That is the point.
psfm(model_name = "GTRE") now reports a persistent scale that has
collapsed to zero, in $sigh_at_bound and $sigr_at_bound and in a
warning naming the surviving scale. Under both estimators: "fiml",
which is the default and packs its results in a separate branch, and
"sml".
GTRE has two persistent components and the likelihood cannot always
separate them in a given sample. When it cannot it merges them: one goes to
zero, the other absorbs its variation and comes back inflated. On simulated
data with both genuinely present, one of the two collapsed in 37% of 87
replications, and sigr -> 0 was three to five times commoner than
sigh -> 0.
Whether firms or periods help depends on the estimator, and the two answers are opposite. Crossed (N, T) designs on the same DGP, regressing the collapse rate on log N and log T:
| estimator | fits | log N | log T |
|---|---|---|---|
| "fiml" (the default) | 600 | -1.66 (0.18) | -0.93 (0.20) |
| "sml" | 251 | -0.68 (0.26) | -1.07 (0.27) |
Under the default, firms buy about 1.8 times what periods do; under
simulated ML it is the other way round. Restricting the "fiml" design to
N <= 200 -- cell for cell the grid the "sml" one used -- leaves it at
-1.64 (0.22), so this is the estimator and not the design; the difference in
the N coefficient is significant (z = 2.8) and the difference in the T
coefficient is not.
So the practical advice under the default is add firms: at T = 4 the
collapse rate falls 68%, 36%, 18%, 6% across N = 50, 100, 200, 400 and
reaches 0% at N = 400, T = 16. Under "sml" that same row is flat (60%,
39%, 41%), which is why a short panel is a dead end there however many firms
it has. ?psfm now gives both.
A user-facing consequence: you can tell in advance. The quantity that
governs whether a panel supports the split is not sigma_h/sigma_r -- that
only decides which scale dies -- but the persistent variance in the firm
means against the transient variance still in them,
[sigr^2 + sigh^2 (1-2/pi)] / ([sigu^2 (1-2/pi) + sigv^2] / T). Across six
designs varying all four scales it orders the collapse rate perfectly
(Spearman -1 under both estimators; logistic coefficient on log(S/N) of
-1.15, SE 0.20): 54% of fits merge at S/N = 0.67, 31% at 1.30, 6% at 11.7.
Below about 1, expect half your fits to merge; it takes roughly S/N > 5 to
get under 10%. ?psfm and the vignette give the formula and a helper.
One limit, stated because it is easy to over-read: those elasticities were
measured at sigma_r = 0.20 against sigma_h = 0.40. On a design where the
random effect dominates instead (0.50 against 0.40, N = 100, T = 6) the
level ordering reverses -- "fiml" collapses on 33% of replications and
"sml" on 20%. Which estimator handles the split better depends on how much
persistent inefficiency there is to find, not on the estimator alone.
This is usually the CORRECT maximum likelihood estimate, not a failure: on one such replication the boundary solution beat the true parameter vector by 3.66 log units. It is the panel counterpart of the cross-sectional wrong-skewness result (Waldman 1982). So it is reported rather than prevented -- bounding either scale away from zero would corrupt exactly the samples where the boundary is the answer.
Your intercept moved with it, and the warning now says so. One rule
covers what a collapse does: it conserves the persistent variance and
relocates the persistent mean into the intercept. Those are two different
destinations -- h is half-normal, carrying both a mean
E[h] = sigh*sqrt(2/pi) and a variance Var(h) = sigh^2 (1 - 2/pi), and
the model has a separate place for each: the intercept takes the mean, since
E[y] = beta0 - E[h] - E[u], and the surviving persistent scale takes the
variance. So the survivor does not come back at sqrt(sigh^2 + sigr^2);
that shortcut predicts 0.447 where the answer is 0.313.
| | surviving scale | intercept |
|---|---|---|
| sigh -> 0 | sqrt(sigr^2 + sigh^2(1-2/pi)) = 0.3133 | -E[h] = -0.3192 |
| sigr -> 0 | sqrt(sigh^2 + sigr^2/(1-2/pi)) = 0.5197 | +0.0955 |
Measured over 1000 replications at each of five sample sizes: 0.2973,
-0.3035, 0.5058, +0.0931 -- within 2.5% to 5.1% with nothing fitted, and all
four slightly below prediction, which is what selection produces rather than
a wrong constant. The practical consequence is that a collapsed fit does not
only mislabel the persistent split; it displaces any level read off coef()
by a known and substantial amount.
And it is wrong skewness, measured between firms. Average the OLS
residuals within each firm -- that average is the persistent part of the
error -- and look at their skew. Over 1000 replications the samples whose
between-firm skew came out the wrong way for a production frontier lost
sigma_h to the boundary 69% of the time against 3.6% for the rest,
with an ROC area of 0.955 for that one statistic, and not one of them
collapsed sigma_r instead. So a collapsed sigma_h is not an accident to
retry with different starting values: the sample carries no between-firm
evidence of persistent inefficiency and the estimate says so.
Read the sign, not a p-value. The firm mean of residuals contains
-ubar_i as well as -h_i, and the mean of T half-normals is skewed the
same way, so the between-firm skew leans negative even at sigma_h = 0 --
82% of the time at T = 4, 63% at T = 10, 53% at T = 30. skewness_test()
is therefore not applicable to between-firm residuals as it stands, and
?psfm says so.
psfm() also gains keep_objective, which retains the likelihood for the
simulated-ML panel models so it can be evaluated away from the optimum. That
is what distinguishes a weakly identified parameter from a badly estimated
one, and it is how the above was established. ?psfm has a worked example.
Robust and clustered standard errors, via sandwich. vcov() returned
the inverse Hessian and nothing else, which is valid only if the likelihood
is correctly specified and the observations are independent. Neither is safe
to assume in applied frontier work, and clustered errors are routine on
firm-level panels.
bread() and estfun() methods are now registered for "sfareg" when
sandwich is installed, which is all that package needs:
fit <- sfm(y ~ x1 + x2, model_name = "NHN", data = d, keep_objective = TRUE)
lmtest::coeftest(fit, vcov. = sandwich::vcovCL(fit, cluster = d$firm))
Two requirements, both of which error clearly rather than returning
something wrong. The fit needs keep_objective = TRUE, since the scores are
differenced from the retained likelihood; and the robust divergence
estimators are excluded, because they do not maximise a log-likelihood and so
have no score. sandwich::vcovHC() does not work and is not expected to --
its corrections are built from hat values, which are undefined for a
nonlinear likelihood whose parameters include variance components.
The scores are central differences with a per-parameter step. Validated two
ways at the optimum: against NHN's hand-derived analytic gradient (agreeing
to 3e-07) and against numDeriv (2.9e-07), with the column sums vanishing
relative to their own scale, which is the identity that makes them scores.
sandwich and lmtest are in Suggests; neither is required to use the
package.
New function skewness_test(): a p-value for wrong skewness. sfm()
already detected the condition -- $wrong_skew, $sigma_u_at_bound, and a
warning -- but could not say whether the skew was wrong by more than sampling
noise. Wrong skew is the single most common reason an applied SFA fit is
meaningless, so that gap mattered.
Two tests: D'Agostino's (1970), which is the one Schmidt and Lin (1984)
appeal to, and Coelli's (1995) M3T. Returns an ordinary "htest", so it
prints like t.test() and carries m3, nobs and wrong_skew alongside
the usual fields. Takes a fitted sfm() object -- using the OLS residuals
the fit recorded, which is what the tests are defined on, rather than the
composed residuals -- or a bare numeric vector.
"agostino" is the default because it holds its size where the asymptotic
form does not. One-sided at a nominal 5%, 4,000 replications of symmetric
residuals:
| n | 25 | 50 | 100 | 400 |
|---|---|---|---|---|
| coelli | 0.031 | 0.037 | 0.043 | 0.050 |
| agostino | 0.045 | 0.046 | 0.046 | 0.051 |
Coelli's is conservative below about n = 100 and the two agree by n = 400.
The D'Agostino implementation agrees with moments::agostino.test to within
1e-10 at n = 30/60/200/1000.
sfm() fits now also carry $ols_residuals, which is what makes the test
exact rather than re-derived from the recorded call.
lcsfm() -- the latent class stochastic frontier (Greene 2005; Orea and
Kumbhakar 2004; Caudill 2003), a new exported entry point alongside sfm(),
psfm(), zsfm() and ttsfm(). J distinct
technologies coexist in one sample, each with its own frontier beta_j and
its own sigma_v,j and sigma_u,j, and which one a firm operates is
unobserved. Class membership follows a multinomial logit with class J as
the reference, constant under "LCM" and parameterized by the second formula
segment under "LCM_Z". The number of classes is n_class, default 2.
Every firm contributes to every class weighted by its class probability, so
this is not the same as fitting separate frontiers to known subgroups -- the
groups are not known, and the uncertainty about them is carried into the
efficiency predictions. post.prob is the n x J matrix of posterior class
probabilities, jlms_class the class-conditional JLMS scores, and jlms
their posterior-weighted average.
It is its own function rather than a model_name inside zsfm(). The two
are related -- zsfm()'s "ZISF" is the restricted two-class case in which
one class has no inefficiency at all -- but a user looking for a latent class
model would not think to look inside a function named for zero inefficiency,
and n_class means nothing to ZISF, whose two components are not free.
Recovery on simulated two-class data, n = 800, true shares 0.60/0.40:
beta_1 = (1.00, 0.48, 0.51) against a true (1.0, 0.5, 0.5),
beta_2 = (2.97, 0.98, 0.21) against (3.0, 1.0, 0.2), scales
0.51/0.98 against 0.5/1.0, fitted shares 0.598/0.402, and 92% of
firms assigned to their true class. On three classes at n = 1200 the
figure is 95.5% after matching labels.
Class labels are identified only up to permutation. Starting values come
from splitting the OLS residuals at their J-quantiles and refitting within
each group, which makes a run reproducible on the same data and avoids the
saddle point that identical starting components would sit at, but it does
not pin labels to any external ordering. Match classes before comparing two
fits. The likelihood-ratio test of J against J+1 is not valid here --
the null puts a class probability on the boundary -- so use AIC()/BIC()
and interpretability instead.
psfm() gave every firm the SAME simulation draws. A single R x 2
Halton block was built once and recycled across all N firms, in three
separate places. sfm() already does the opposite and says why: Train (2002,
p. 228) attributes Halton's advantage to its coverage and to the negative
correlation it induces across observations, and the second half only exists
when observations get different blocks. Sharing one block also makes every
firm's simulation error the same realization, so it cannot average out as N
grows. Firm i now gets its own contiguous block.
Measured, GTRE by simulated ML, 30 replications at each of N = 50/100/200,
old and new on identical seeds so the comparison is paired: sigh improves
on 62 of 87 pairs (Wilcoxon p = 0.015) and sigr on 48 of 87
(p = 0.048), while beta1 -- the control, since the frontier slope is not
what the draws integrate over -- is unaffected at 46 of 87 (p = 0.41).
Collapses of sigh to exactly zero halved, 11/87 to 6/87.
The improvement grows with N, which is the point: 16/28 at N = 50
(p = 0.58), 21/29 at N = 100 (p = 0.25), 25/30 at N = 200 (p = 0.012). That
is the shape the mechanism predicts.
It does not fix the GTRE convergence failure. RMSE(N=50)/RMSE(N=200) for
sigh moves from 0.94 to 1.08 where root-n requires 2.00. Shared draws were
a contributor, not the cause.
rand.gtre randomized the draws in a way that destroyed what they were
for. The old code drew 9999 random permutations of Halton dimension 1 and
kept whichever correlated least with dimension 2. A two-dimensional Halton
sequence is worth having because the pairs cover the unit square evenly;
permuting one column preserves each margin, randomizes the pairing, and
leaves joint coverage no better than random. Measured at R = 150 with
primes 2 and 3 and 1000 discarded: joint discrepancy 0.0217 before the
shuffle and 0.0550 after, against 0.1400 for purely random pairing -- to
remove a correlation of 0.008 that was never a problem.
rand.gtre now applies a uniform shift modulo 1 (Tuffin 1996; Train 9.3.4),
which moves the lattice without disturbing its structure: same measurement,
0.0250. Results change for anyone who passed rand.gtre.
psfm(model_name = "KSS") -- Kneip, Sickles and Song (2012). The firm
effect is a smooth function of time lying in an L-dimensional space whose
basis is estimated from the data:
alpha_it = sum_r theta_ir * g_r(t). Nothing on CRAN implements this, so it
is a differentiator rather than catch-up.
It makes explicit what the other time-varying estimators are special cases
of: SSFE is L = 1 with a constant basis, LS is L = 1 with the basis
free, CSS is L = 3 with the basis fixed to {1, t, t^2}, and KSS
estimates both. Follows the original's three steps -- cross-sectional
centering by period, then smoothing of each firm's residual trajectory
followed by an eigendecomposition of their empirical covariance, then
loadings by least squares -- iterated with beta. Balanced panels only,
which is what the estimator is defined on; it refuses rather than quietly
fitting something else.
Recovers the true rank exactly at L = 1, 2, 3 on simulated data
(n = 110, T = 10), with the estimated basis spanning the true factor space
and sigma_v recovered to within 3%.
Both tuning choices are exposed. kss_smooth defaults to GCV; kss_L
defaults to the Bai and Ng (2002) IC_p2 criterion, which is not the
original paper's threshold rule -- said plainly in ?psfm because a
different rule can select a different L. $kss$eigenvalues is returned so
the choice can be inspected.
The automatic search is capped at floor(T/2), and warns if it selects that
cap. IC_p2 works because the residual variance flattens out once the real
factors are in; on a short panel, letting L reach T - 1 lets the factors
span nearly the whole time dimension, so the variance collapses and no
penalty competes. Measured at n = 100 over five seeds and two true ranks, an
uncapped criterion returned its maximum on all ten designs at T = 6 and
again at T = 8. With the cap it recovers the true rank exactly from T = 8
upward. An explicit kss_L is the user's own call and is not subject to the
cap.
One trap worth knowing: KSS's period centering means its alpha_hat
carries no period mean while CSS's and LS's do, so the two differ by a
per-period constant by construction. Correlating them directly understates
the agreement badly (0.93 against 0.99 on a design where they agree). Compare
u_hat, which is a within-period contrast and invariant to the shift.
psfm(model_name = "SSRE") and "SSCRE" complete the Schmidt-Sickles
(1984) family that SSFE starts: the GLS estimator, and its
correlated-random-effects correction. SSRE is more efficient than SSFE
and identifies time-invariant regressors, but assumes the effects are
uncorrelated with the regressors -- in a frontier model, that a firm's
inefficiency is unrelated to its input choices.
SSCRE adds the within-firm means of the time-varying regressors (Mundlak
1978), which models that correlation rather than assuming it away. Its
slopes then equal SSFE's within slopes exactly, which the tests pin at
1e-6; the coefficients on the added means (reported with a .mean_ prefix)
are the Mundlak form of the Hausman test.
Two classical time-varying panel estimators: psfm(model_name = "CSS")
and psfm(model_name = "LS"). Cornwell, Schmidt and Sickles (1990) and
Lee and Schmidt (1993). Both sit between "SSFE", which holds inefficiency
fixed over time, and the ML panel models, which let it move only by assuming
a distribution for u and a decay path for it: these let the firm effect
itself vary with time and read inefficiency off it, assuming no distribution
at all.
They differ in what they spend on flexibility. CSS gives every firm its
own quadratic in time (3N parameters; firms may overtake one another;
T_i >= 4 required before a firm contributes to beta). LS imposes one
common temporal pattern scaled per firm, alpha_it = delta_t * alpha_i
(N + T - 1 parameters; the ranking of firms cannot change, only the
spread). Cross-over versus no cross-over is the substantive choice between
them.
CSS residualizes on each firm's own (1, t, t^2) rather than building
3N dummy columns, and refuses a regressor spanned by those quadratics (a
pure time trend, say) instead of reporting it as a small coefficient. LS
is a rank-one factor model and is fitted by alternating least squares, so an
unbalanced panel needs no special-casing. On simulated data at N = 90,
T = 8, LS recovers the temporal pattern with correlation > 0.99 to truth
and the firm scales at 0.995; both recover the frontier slopes and the noise
standard deviation.
Neither is maximum likelihood, so neither carries an $opt component and
logLik()/AIC()/BIC() return NA with a warning, as for SSFE.
The wrong-skew boundary report fired on a knife-edge. sfm() flags
sigma_u_at_bound when the one-sided scale collapses relative to the
residual SD, and the threshold was 1e-3. On the design used to test it the
collapsed sample sits at 8.6e-4 and the interior ones at 0.47-0.71, so
the threshold cleared the case that matters by a factor of 1.16 and the
others by 470. An optimizer's stopping point moves by more than 16% between
BLAS implementations: the same fit reported TRUE locally and FALSE on
another platform. Now 1e-2, which sits in the middle of a gap spanning
nearly three orders of magnitude and leaves every sample at least 11x clear.
This only ever changes the SECOND of two conditions -- the warning also
requires wrong skew (m3 >= 0), which is a deterministic function of the
data -- so the practical effect is that the report is stable across
platforms rather than that it fires more often.
nobs() counted rows SUPPLIED rather than rows USED, so BIC() was
computed against the wrong n whenever any row was dropped for
missingness. nobs.sfareg() re-evaluated the data argument of the
recorded call, which of course returns the full frame; data_proc2() has
always dropped incomplete cases before fitting. On a 300-row frame with 7
missing values nobs() returned 300 against 293 actually used. The fitted
object now records the effective count directly, and nobs() falls back to
the length of a stored per-observation vector before it resorts to
re-evaluating the call.
sfm(model_name = "NE") could return a positive log-likelihood with
sigma_u = 0 and a divergent sigma_v. The log-density was built as
-log(sigma_u) + log Phi(z) + eps/sigma_u + sigma_v^2/(2 sigma_u^2). Using
pnorm(log.p = TRUE) for the middle term is correct as far as it goes, but
that term and the tilt both diverge like z^2/2 with opposite signs as
sigma_u -> 0, and their sum is a catastrophic cancellation. At
sigma_v = 587, sigma_u = 1e-7 the two are -/+1.7252e19, where
consecutive doubles are 2048 apart, so the sum came back as rounding noise --
positive noise, which the optimizer then maximized by running sigma_u to
its lower bound.
This was not rare and it was silent: scanning 150 samples at
lambda = 0.75, N = 100, one returned sigma_v = 587.4, sigma_u = 0,
logLik = +468992, and another sigma_v = 3.4e15, logLik = +7.9e30, both
as ordinary sfareg objects with no error and no warning. Across a
12-cell design at 1,500 replications the rate was 0.74%, reaching 4.4% at
lambda = 0.5, N = 100.
New .log_phi_tilt() (matrix_utils.R) does the cancellation analytically
instead. Because
z^2/2 = eps^2/(2 sigma_v^2) + eps/sigma_u + sigma_v^2/(2 sigma_u^2), the
tail expansion of log Phi(z) cancels the tilt exactly and leaves only
-eps^2/(2 sigma_v^2) - log(2 pi)/2 - log(-z) + log1p(-1/z^2 + 3/z^4 - ...),
in which no large intermediate is ever formed. It agrees with the previous
expression to 8.5e-13 over 75,000 evaluations wherever the previous one was
trustworthy, integrates to 1, and takes this failure rate to 0 in 2,200
fits.
To be precise about what that last number does and does not cover: it is the
rate of the catastrophic-cancellation failure -- a positive log-likelihood
with sigma_u at its bound -- and nothing else. It is not a claim that NE
fits are uniformly well behaved. convergence/CONVERGENCE_STATUS.md records
a separate and still-open defect at the other end of the parameter space: at
lambda = 3.33 the root-n sweep shows a tail of 2-4% of replications landing
far from the truth on some parameter, and unlike the rest it does not shrink
with n. There the scales are near the truth and the frontier SLOPES break,
with beta1 and beta2 trading against each other, so it is a different
failure from this one and is not addressed by this fix. The two live in
different regions -- cancellation bites at LOW lambda, the tail at high -- and
the median fit is unaffected in both.
"NGE" carried the identical defect and is fixed the same way. Its
likelihood is a difference of two exponentially tilted Gaussians, both with
the same structure as "NE"; both now go through .log_phi_tilt(), in the
likelihood and in the post-estimation efficiency block. Agreement with the
previous expression is 7.1e-10. "TSL", "TTNE" and "NE_Z" share the
structure and have not been checked.
"NE" starting values now come from a bias-corrected moment estimator
(new R/ne_start.R) rather than the flat sigma_u = sigma_v = 0.1 that
start_cs() hands every cross-sectional model. Minus the mean negative OLS
residual has an exact asymptotic bias factor h(lambda); .ne_start()
divides it out, shrunk by the share of MSE the bias accounts for. Over an
88-cell design at 2,000 replications, total MSE of the start against the
truth: 5.49 for the uncorrected version, 2.86 for the third-moment (COLS)
inversion, 1.88 for their minimum, 1.82 for this one. Effect on the fitted
MLE is smaller and confined to lambda >= 1.5, where it lowers MSE of
sigma_u-hat by 5-51%; below lambda = 1 it changes nothing. Derived in
Bernstein, Parmeter and Wright, "Starting Values for the Normal-Exponential
Stochastic Frontier Model".
Note this leaves the flat 0.1 start in place for every other
cross-sectional model, which PROJECT_STATUS.md still lists as open.
sfm() now reports when sigma_u sits on the zero boundary under wrong
skewness, instead of saying nothing. A one-sided scale at zero when the OLS
residuals are skewed the wrong way is the correct maximum likelihood
estimate -- the Type I failure of Olson, Schmidt and Waldman (1980) -- not a
numerical problem. The estimator = "cols" path has warned about this for
some time; the likelihood path did not, so a user got a boundary fit with no
explanation.
New $wrong_skew, $sigma_u_at_bound and $residual_m3 components, plus a
warning that says what the boundary means and that the efficiency scores are
uninformative there. Measured at lambda = 0.75, N = 100 over 600 fits:
13.2% of samples are wrongly skewed, 17.7% of those put sigma_u on the
boundary, and not one correctly skewed sample does. That is why the fix
is to report the boundary rather than bound sigma_u away from zero -- a
bound would corrupt precisely the samples where the boundary is the answer.
"NE" no longer emits a stream of NaNs produced warnings. Its
likelihood had no guard against a non-positive scale, so the optimizer
probing sigma_u <= 0 evaluated log(sigma_u) and warned -- 17 times in a
single fit on a wrongly skewed sample, burying the warning that mattered. It
now returns a large finite penalty for out-of-domain parameters, as "NLN"
and "NW" already do. No estimate changes: this bounds the objective's
domain, not the estimate.
"NGE" aborted outright on about 7% of small samples, and "NU" could
have. Both guarded their scale parameters with
return(.Machine$double.xmax). optim() differences the objective to form a
gradient, and differencing 1.8e308 overflows to a non-finite value, so the
fit died with non-finite finite-difference value instead of the optimizer
being steered away. Measured before the change: 3 of 45 "NGE" fits at
N = 150 failed this way, including at sigma_u = 1, sigma_v = 0.3. Both now
use the same finite penalty the other branches use -- 0 of 45 failures after,
and "NU" estimates are bit-identical across 45 fits, since the penalty only
ever applies where the parameters are already inadmissible.
ttsfm() gains z_link too, completing the cross-sectional half of the
variance-determinant inconsistency. Two-tier models placed z'delta on the
standard deviation, like sfm() and unlike psfm() and the competing
packages; z_link = "var" aligns them. The default is "sd", and fits under
it are byte-identical to before -- checked against the previous build with
the particle swarm seeded, for both "TTNE" and "TTHN".
Both tiers move together: delta on zu and on zw each halve exactly
between the links, while the maximised log-likelihood and the frontier
coefficients are unchanged.
marginal_effects() now works on panel fits, psfm(model_name = "TRE_Z")
and psfm(model_name = "GTRE_Z"). These place z'delta on the variance,
and the link is read from the fit rather than assumed, so their effects are
directly comparable with an sfm() fit -- including one fitted with the new
z_link = "var".
"GTRE_Z" separates persistent inefficiency from transient, and both
are reported: component = "u" for the transient block and
component = "h" for the persistent one. Each block is located by name
rather than by position, because the sigma_h coefficients are the
trailing ones -- a positional rule would have silently reported sigma_h
effects under a sigma_u label. Output columns carry the component
(dE_u.dz against dE_h.dzp), so a table cannot be misread once separated
from the call that produced it. The default is "u", so nothing existing
changes.
sfm() gains z_link, which fixes a real trap in comparing _Z fits
across the package. sfm()'s "NHN_Z"/"NE_Z" put the
variance-determinant linear predictor on the standard deviation,
sigma_u = exp(z'delta), while psfm()'s "TRE_Z"/"GTRE_Z" put it on the
variance, sigma_u = sqrt(exp(z'delta)) -- as do the competing packages.
Since exp(eta) = exp(eta/2)^2 the two fit the same model and return the same
sigma_u, the same log-likelihood and the same marginal effects, but delta
under the SD link is exactly half delta under the variance link. Reading
a delta from one family as if it came from the other therefore doubles or
halves every reported effect.
z_link = "var" puts an sfm() fit on the same footing as psfm() and
sfaR. The default is "sd", so no existing result changes. The
efficiency predictor and the z_spec that marginal_effects() reads both
follow whichever link was used.
sfm()'s simulated-ML models "NLN" and "NW" now estimate the composed
density from two proposals at once, which fixes an accuracy defect that no
single-proposal scheme can. The integrand is a product of a normal kernel of
width sigma_v centred at u = -eps and the inefficiency density. Drawing
from the inefficiency quantile -- what both models did through 1.1.5 -- misses
the kernel whenever it is the narrow factor: at eps = -4.7 with
sigma_v = 0.3 the spike needs a draw beyond the 0.9999 quantile, while 400
draws reach only about 0.9975, and five observations in 800 carried 72% of
"NW"'s total error. Drawing from the noise instead cures that case and
breaks the mirror one, wherever the inefficiency density is the narrow factor.
Both draws are now taken -- half the count each -- and combined by the balance heuristic of Veach and Guibas (1995), each draw weighted by which proposal was likelier to have produced it. There is no selection rule to get wrong. Total error against adaptive quadrature, over an 80-cell parameter grid on a 300-observation sample, restricted to the cells an optimizer can reach:
| cells within ... of the optimum | 1.1.5 scheme | now |
|---|---|---|
| "NW", 200 | 220, worst cell 172 | 16, worst cell 5.7 |
| "NW", 500 | 3633, worst cell 1910 | 40, worst cell 11 |
| "NW", 1000 | 17900, worst cell 8553 | 128, worst cell 64 |
| "NLN", 200 | 233, worst cell 120 | 11, worst cell 1.0 |
| "NLN", 500 | 3564, worst cell 1553 | 49, worst cell 8.8 |
| "NLN", 1000 | 15962, worst cell 6639 | 205, worst cell 63 |
(Each scheme is measured at its own draw rule, so "NW"'s 1.1.5 column uses
the 400 draws that version would have taken and the new column 200.)
"NW"'s draw rule drops from max(400, 8*sqrt(n)) to the max(200, 3*sqrt(n)) "NLN" already used, so at n = 1500 an "NW" fit costs about
what it did in 1.1.5 (roughly 50 seconds) while being far more accurate;
"NLN" costs about 1.5x more, since its draw count did not fall. Nsim
still means the total number of draws.
The two models now also converge at the same rate, so the old advice that
"NLN" needs far more draws than "NW" no longer applies -- that was a
property of the integration scheme, not the model. At n = 3000, measured
against quadrature at the true parameters, total simulation error is 3.4
log-likelihood units for "NLN" and 3.0 for "NW" at the default, falling
to 0.07 and 0.06 at Nsim = 6400. The comparable figure for the old scheme
was 226.8.
The accuracy gain shows up in the estimates. Over 12 replications at n = 800, bias and RMSE against the true parameters:
| | sigma_v | sigma_u | shape | x1 |
|---|---|---|---|---|
| "NW" 1.1.5 | +0.036 / 0.066 | -0.149 / 0.237 | -0.160 / 0.293 | -0.021 / 0.043 |
| "NW" now | -0.006 / 0.044 | -0.024 / 0.112 | +0.001 / 0.111 | +0.001 / 0.045 |
| "NLN" 1.1.5| -0.012 / 0.022 | +0.005 / 0.101 | +0.012 / 0.166 | +0.003 / 0.037 |
| "NLN" now | -0.004 / 0.035 | +0.003 / 0.080 | +0.005 / 0.131 | +0.002 / 0.038 |
"NW" is the clear case: its sigma_u and shape bias essentially vanish and
their RMSE roughly halves. "NLN" improves more modestly, and its sigma_v
RMSE is slightly worse -- at 12 replications that difference is within noise
and should not be read either way.
The efficiency predictor reuses the weights the likelihood computed, so
u_hat cannot drift away from the density that was maximised.
sfm() gains user control over its simulated-ML draws, following the
practice in Train (2002, ch. 9) and matching the options sfaR exposes:
sim_type ("halton", "sobol", "torus", "uniform"), antithetics,
sim_burn, sim_scrambling, sim_prime and sim_seed. These affect the
models fitted by simulation, "NLN" and "NW".
"sobol" with sim_scrambling 1-3 is the scrambled sequence of Bhat
(2003), which removes cross-dimension correlation while keeping the
coverage. "uniform" exists as a baseline to measure against, not as a
recommendation -- Bhat (2001), quoted by Train, found 100 Halton draws more
precise than 1000 pseudorandom ones. antithetics = TRUE takes half the
draws from the sequence and creates the rest as mirror images (Hammersley
and Morton 1956), which costs nothing.
Defaults reproduce the previous behaviour exactly -- Halton, no antithetics, 1000 leading elements discarded -- and this is asserted by a test that rebuilds the old construction and compares byte for byte, so no existing result moves.
The draw construction is now one shared internal function,
.sml_draws(), rather than being written inline at each site. That matters
beyond tidiness: the sites had drifted apart, with sfm() giving each
observation its own contiguous block of the sequence while psfm() hands
every firm the same one. Train (ch. 9) attributes the value of a
low-discrepancy sequence to its coverage and to the negative correlation
it induces across observations, and the second only exists when units get
different blocks. The shared constructor always blocks by unit and takes a
dim argument for the multi-dimensional panel case, so the other entry
points can adopt it directly.
New model sfm(model_name = "TSL"): the normal / truncated skew-Laplace
frontier (Wang 2012). Inefficiency has the signed-mixture density
f(u) = ((1+lambda)/(sigma_u(2*lambda+1))) [2 exp(-u/sigma_u) - exp(-(1+lambda) u/sigma_u)] on u >= 0 -- the second exponential enters
negatively -- so it nests the exponential model as lambda -> 0 while
allowing a non-monotonic inefficiency density. Reports
(sigv, sigu, lambda) plus the frontier, with both u_hat and
exp_u_hat.
The composed density is evaluated as a difference in log space rather than
by subtracting the raw terms. Both carry sigma_v^2/(2 sigma_u^2), which
overflows once sigma_u is small relative to sigma_v -- on a routine grid
spanning sigma_u in {0.05, 0.1, 0.3} and sigma_v in {0.3, 1} the direct
form returns NaN or Inf at 9 of 30 points, including sigma_u = 0.05,
sigma_v = 1, eps = 0, which is not an extreme point. The log form is
finite throughout.
Verified two ways: against numerical convolution of the implied u density
with the normal noise (agreeing to 8 decimals across lambda in
{0.5, 1.5, 4}), and against an independent implementation, which it matches
to an identical log-likelihood and identical coefficients on 8 of 8 samples.
lambda is a shape parameter and is the least sharply identified of the
three -- read its t-value before interpreting it.
data_gen_cs() gains lam_tsl and the matching u_tsl / y_pcs_tsl
columns. The signed mixture cannot be drawn by picking a component, so it
is drawn by rejection off its first exponential; the acceptance probability
never falls below 1/2.
New function efficiency_ci(): Horrace and Schmidt (1996) confidence
intervals for individual inefficiency. sfa has always reported point
predictions of u_i -- u_hat (Jondrow et al. 1982) and exp_u_hat
(Battese and Coelli 1988) -- and nothing about how sharply either is pinned
down. Both are posterior means of u given the composed residual, and the
posterior is a truncated normal in closed form, so the interval costs no
estimation beyond the fit itself. Returns a data frame of
u_lower/u_hat/u_upper and the corresponding technical-efficiency
bounds, at any level.
Available for sfm()'s "NHN", "NHN_Z", "NE" and "NTN", whose
posterior really is a truncated normal. Other models have a posterior of a
different shape, and the function says so by name rather than returning a
number the formulas do not support.
These intervals condition on the fitted parameters: they do not narrow as
n grows, because they measure the irreducible difficulty of splitting one
residual into noise and inefficiency, not estimation uncertainty. That is
worth knowing before reading an efficiency ranking closely -- on a routine
200-observation half-normal fit the median 90% interval for a single unit
spans about 0.49 in efficiency, against a total spread of 0.85 across all
the point predictions.
npsfm(method = "SZ") now solves its DEA step in the package. The
output-oriented envelopment program is computed directly, as one linear
program per unit, by the new internal .dea_out(), covering all four
returns-to-scale settings. Suggests gains lpSolve, a general-purpose
linear-programming solver (pure C, no system libraries).
Results are unchanged, and were verified against an independent reference implementation over all four returns-to-scale settings at n = 30, 80 and 150, with one to three inputs and one to two outputs: agreement to 3e-12 or better throughout.
New function marginal_effects(): what the _Z models' delta
coefficients actually mean for inefficiency. For a fit whose inefficiency
scale depends on covariates, it returns the per-observation
d E[u]/d z_k and d Var[u]/d z_k, with the average marginal effects
attached. These are what applied papers report; a delta on its own is a
coefficient in a log link for a scale parameter and is not interpretable in
the units of either u or z.
This also defuses a long-standing trap. sfm()'s "NHN_Z"/"NE_Z" put the
linear predictor on the standard deviation and psfm()'s "TRE_Z"/
"GTRE_Z" put it on the variance, so the same delta means different
things in the two families -- the marginal effects differ by a factor of
two. The effect is on the scale of u either way, so reporting it rather
than the coefficient removes the ambiguity. The returned object records
which convention was used.
Currently covers sfm()'s "NHN_Z" and "NE_Z"; psfm()'s panel _Z
models are not wired up yet. Standard errors are deliberately not reported --
see ?marginal_effects for why.
"NHN_Z" and "NE_Z" fits now carry a $z_spec component (the
variance-determinant design, its coefficients, and the link convention),
which is what marginal_effects() reads.
sfm(model_name = "NNAK") now starts from the method of moments. The
normal-Nakagami likelihood has sigma_u -> 0 as a genuine attractor, and
the previous hard-coded start of sigma_u = sigma_v = 0.1 sat next to it.
Over twelve samples at n = 3000 with a true sigma_u of 1, the old start
drove sigma_u to 0.0013 and to 0.0000 on two of them -- inefficiency
vanishing altogether, with the frontier intercept pushed negative -- for a
log-likelihood 44 and 45 points worse than the new start reaches. The new
start was never worse on any of the twelve and was strictly better on seven,
for a mean gain of 7.9 log-likelihood points.
The construction follows FronPy (Stead 2024, Journal of Productivity
Analysis, the paper these closed forms come from): invert the half-normal
moment equations for the two scales and shift the frontier intercept up by
the implied E[u]. The half-normal is the right auxiliary because it is the
m = 1/2 member of the Nakagami family. The shape still starts at 0.5, as
before. When the residuals are skewed the wrong way and the moment equations
have no admissible solution, the old constants are used as before.
sfm(estimator = "mols") is accepted as a synonym for
estimator = "cols". What the package computes under that name is the
modified OLS moment estimator of Olson, Schmidt and Waldman (1980) -- it
inverts the second and third central moments of the OLS residuals -- and not
Winsten's corrected OLS, which shifts the intercept by the largest
residual and estimates no variance parameters. The documentation now names
it MOLS and draws the distinction. "cols" is unchanged and not deprecated.
Fitted "sfareg" objects from the four models above carry a new
$u_posterior component (mu_star, sigma_star), which is what
efficiency_ci() reads.
The intro_to_psfm vignette now builds in a fraction of the time. At
CRAN's request, the models it fits are sized as toy illustrations rather
than as estimation exercises: the simulated panel is 70 firms over 6
periods instead of 100 over 10, the simulated-ML fits draw 50 Halton
points via halton_num rather than the default
ceiling(sqrt(nrow(data))) + 100, and the psfm_bootstrap() example
uses 5 replications on a 30-firm panel instead of 10 on 60. Vignette
rebuild time falls by roughly a factor of five. The reported estimates
therefore differ from previous versions, and the vignette now says
plainly that they are not to be read as a serious fit.
Every model fit in that vignette now sets rand.gtre and rand.psoptim.
With PSopt = TRUE the particle-swarm stage draws from the session's RNG,
so the vignette's results previously changed from one build to the next;
they are now reproducible.
No change to any R code, to NAMESPACE, or to the documented interface.
psfm(model_name = "GTRE") now defaults to full information maximum
likelihood. The four ways of fitting the four-component GTRE model were
four separate model_name values, which made them look like four different
models rather than four routes to the same one. They are now selected with an
estimator argument, in the same spirit as sfm()'s estimator = c("mle", "cols"):
"fiml" (the default) -- full information ML through the closed-skew-normal
representation. Deterministic; requires a balanced panel."sml" -- simulated ML over Halton draws. Handles unbalanced panels. This
is what "GTRE" meant through 1.1.3."seq1", "seq2" -- the two-step moment estimators.Scripts that pass model_name = "GTRE" therefore get a different estimator
than they did, and are warned once per call. Pass estimator explicitly to
silence it. The names "GTRE_FML", "GTRE_SEQ1" and "GTRE_SEQ2" are
unchanged and still select the same three routes directly.
On an unbalanced panel "fiml" cannot be fitted. Taking the default
warns and falls back to "sml", because erroring would make "GTRE"
unusable by default on a whole class of data; asking for "fiml" explicitly
errors instead of silently fitting something else.
psfm(model_name = "PL80_MVTN") -- Pitt and Lee's (1981) Model III, the
multivariate truncated normal panel likelihood from their Appendix 2.
Inefficiency varies over time and is correlated within a firm:
u_i = (u_i1,...,u_iT)' is drawn from a T-variate normal truncated to the
negative orthant, where the existing "PL80" holds inefficiency fixed over
time.
Pitt and Lee derived this likelihood and then set it aside, writing that it
"is difficult to evaluate since the quantities P0 and P(y_i - x_i beta)
involve T-dimensional numerical integrals", and estimating Model III by
Zellner SUR instead. Those quantities are multivariate normal orthant
probabilities; mnormt::sadmvn() evaluates one in about 3 ms at T = 6, so a
likelihood evaluation costs roughly N+1 of them -- about 0.3 s at N = 100,
and about 12 s for a whole fit at N = 80, T = 4. What was intractable in
1981 is merely slow now.
Sigma is parameterized as equicorrelated, sigma_u^2[(1-rho)I + rho 11'],
costing two parameters. Unrestricted Sigma costs T(T+1)/2 -- 21 at T = 6,
55 at T = 10 -- every one identified only through orthant probabilities. The
equicorrelated form captures what the general Sigma was for: dependence of a
firm's inefficiency across periods, with rho = 0 giving independence over
time and rho -> 1 approaching the time-invariant "PL80". Because the form
is equicorrelated, every matrix quantity is closed form (Sherman-Morrison,
verified against brute force to 1e-17), so only the orthant probabilities are
numerical.
Requires a balanced panel with T >= 2; an unbalanced one errors and points
at "PL80". data_gen_p() gains y_pl_mvtn and u_mvtn columns and a
rho_mvtn argument to test it, generated last so the RNG stream feeding every
existing column is untouched. Registered in the convergence framework as
PL80_MVTN. See ?PL80_MVTN.
Validated before wiring: profiling the likelihood one parameter at a time
puts the minimum of the negative log-likelihood at the truth for all five
parameters, and maximizing it recovers (sigma_v, sigma_u, rho) =
(0.32, 0.78, 0.45) against a truth of (0.30, 0.80, 0.50) at N = 300.
npsfm(), nonparametric stochastic frontier models. A fifth entry point,
alongside sfm(), psfm(), zsfm() and ttsfm(), for frontiers whose shape
is estimated by kernel regression rather than assumed linear. Two estimators:
method = "FLW" -- Fan, Li and Weersink (1996). Fits E[y|x]
nonparametrically, then recovers the scale parameters from the residuals,
by maximizing their concentrated likelihood in lambda for
dist = "hn" and by inverting central moments for dist = "exp",
"gamma" and "unif".method = "SVKZ" -- Simar, Van Keilegom and Zelenyuk (2017). Three
local-linear regressions give sigma_u(x) and sigma_v(x) pointwise, so
both variance components vary with the covariates. No optimizer runs.
Normal-half normal only.method = "PSZ" (alias "KPST") -- Park, Simar and Zelenyuk. Local
maximum likelihood: the frontier and both log variance components get
local-linear expansions and the kernel-weighted normal-half normal
likelihood is maximized in those 3(k+1) parameters, once per observation.method = "MY" -- Martins-Filho and Yao. Iterative local likelihood,
alternating local frontier fits with a global update of (lambda, sigma).method = "SZ" -- Simar and Zelenyuk (2011). Passes an existing smooth
frontier through an output-oriented DEA to impose monotonicity and
convexity."PSZ" and "MY" run one numerical optimization per observation (for
"MY", per observation per iteration), so they are one to two orders of
magnitude slower than "FLW". Both are seeded from an "FLW" fit.
Ported from Christopher Parmeter's research scripts. Results return as class
"npsfareg" rather than "sfareg": there is no parameter vector with
standard errors, so coef(), vcov() and logLik() would have nothing to
return. fitted(), residuals(), nobs(), print() and summary() are
provided.
Kernel regression comes from the np package, added to Suggests rather
than Imports -- nothing else in sfa needs it, and npsfm() checks for it
and stops with an install instruction if it is absent.
A correction worth recording, because it is easy to repeat: the two
local-likelihood estimators maximize the composed-error likelihood, in
which the local intercept is the frontier m(x) itself. They therefore take
no half-normal mean shift, unlike the least-squares methods, whose first
stage estimates E[y|x] = m(x) - E[u] and does need one. Applying the shift
to "PSZ"/"MY" biases the whole frontier up by about E[u]; mean absolute
frontier error at n = 300 fell from 0.372 to 0.116 ("PSZ") and 0.410 to
0.070 ("MY") once it was removed.
Against a simulated nonlinear frontier with sigma_u = 0.6, sigma_v = 0.25,
both least-squares estimators converge as n grows (6 replications at each size):
FLW recovers sigma_u = 0.560, 0.539, 0.595 at n = 150, 300, 600, and
SVKZ 0.477, 0.506, 0.559, with mean absolute frontier error falling from
0.102 to 0.043 and 0.206 to 0.079 respectively. SVKZ's downward bias at
small n is the wrong-skew floor: the share of observations whose local
third moment has the wrong sign, and whose sigma_u(x) is therefore set to
zero, falls from 21.6% to 0.3% over that range.
sfa_diagnostics(), plot() for "sfareg", and convergence reporting.
The numerical hardening was already in place -- staged minimizer, clipping
constants, analytic gradients where they exist -- but nothing was reported
back. A fit carried optim()'s convergence code, message, evaluation counts
and Hessian, and print()/summary() showed none of it, so a fit that
stopped on the iteration cap printed exactly like a converged one.
sfa_diagnostics() returns the convergence code and what it means, the
eigenvalue spectrum and condition number of the Hessian, whether it is
positive definite, which parameters load on its flattest direction, the
implied parameter correlations, and -- with keep_objective = TRUE -- the
gradient at the reported optimum. plot() draws the Hessian spectrum, the
correlation matrix, a likelihood slice per parameter, and the gradient.
print() and summary() now report the convergence code.
The code by itself is not diagnostic and is not treated as though it were.
Across NHN, NE and NTN at n = 150, 500 and 1500, code 52
("ABNORMAL_TERMINATION_IN_LNSRCH") appears routinely alongside a maximum
relative gradient of ~1e-6 and a positive definite Hessian: the staged
minimizer had already converged and L-BFGS-B could not step away from the
optimum. The same code on NTN at n = 150 came with a relative gradient of
5e+07 and an indefinite Hessian, a real failure. The verdict therefore
combines the code with the gradient and the Hessian, and distinguishes
benign from unverified (a line-search code with no objective retained, so
no gradient to settle it) from failure. Code 1, the iteration limit, is
never treated as benign.
On a single NNAK fit the report reproduces what the convergence sweeps
found only across replications: mu and sigu correlated at 0.998 -- the
documented ridge -- with the flattest Hessian direction loading on exactly
that pair.
When the Hessian is singular enough that some parameter has no usable variance, the correlation report drops those parameters and names them, rather than vanishing as a whole -- a diagnostic for ill-conditioning should be most informative exactly when conditioning is worst, not least.
sfm(keep_objective = TRUE) stores the likelihood on the fitted object so
the gradient and likelihood slices can be computed after the fact. Off by
default: a closure carries its enclosing environment, so a fit saved with one
serializes the estimation data too (about 38 KB to 1.7 MB on a 200-observation
example). Everything else sfa_diagnostics() reports works without it.
Model names are matched without regard to case. match.arg() is case
sensitive, so psfm(model_name = "gtre") used to fail with a list of valid
names that visibly contained what the user had just typed. All five entry
points now fold case, for model_name and for npsfm()'s method. No entry
point has a case collision among its choices -- sfm()'s "THT" and "tHN"
differ in more than case -- so the canonical spelling is always recoverable.
Exact matches beat partial ones, which matters because "GTRE" is a prefix
of four other names and must resolve to itself. Genuinely ambiguous partials
("GTRE_S", between GTRE_SEQ1 and GTRE_SEQ2) are still rejected rather
than guessed at, and an unrecognized name now suggests the two closest valid
choices instead of listing everything.
psfm()'s optimizer iteration defaults have been raised, from
maxit.bobyqa = 100, maxit.psoptim = 10, maxit.optim = 10 to
5000, 100 and 1000. maxit.nlminb is now an argument (default 500);
it was previously hard-coded at 200 in the "GTRE_FML" branch and 500
elsewhere, and could not be set from the call.
The old caps were binding rather than merely economical. The
"K1990"/"K1990modified" code already carried a note that 100 bobyqa
evaluations left its seven-parameter fits several log-likelihood units short
of the optimum purely on the iteration cap, and "GTRE_FML" at N = 500,
T = 10 roughly halves its root-mean-square error against known true values
once the caps are lifted (0.0080 to 0.0037 and 0.0179 to 0.0073 on two
draws), for about 1.5 times the run time.
Existing scripts will get more accurate estimates and slower fits. Pass the old values explicitly to restore the previous behaviour.
psfm(model_name = "TFE") now fits a different estimator. Through
version 1.1.2 the name "TFE" selected Chen, Schmidt and Wang's (2014)
within maximum-likelihood estimator. It now selects Greene's (2005) true
fixed effects estimator, which is what the name means in the literature.
The Chen-Schmidt-Wang estimator is unchanged and is now
model_name = "TFE_WMLE".
Scripts written against 1.1.2 or earlier that pass "TFE" will silently get
a different estimator, so psfm() issues a warning whenever "TFE" is used.
To reproduce earlier results, change the name to "TFE_WMLE".
Corrected ordinary least squares, sfm(estimator = "cols"). The moment
estimator of Olson, Schmidt and Waldman (1980). OLS is consistent for the
slopes of a composed-error frontier whatever the one-sided distribution --
only the intercept is biased, by E[u] -- so COLS keeps the OLS slopes,
inverts the central moments of the OLS residuals for the scale parameters,
and shifts the intercept up by the implied E[u].
Implemented for "NHN", "NE" and "NG"; other models error, because the
moment inversion is distribution-specific. No optimizer runs and the result
is deterministic, which makes it a natural robustness check against a
maximum-likelihood fit that may have settled at a local optimum.
Wrong-skew samples are reported rather than absorbed: a production frontier
implies a negative third central moment, and when a sample comes out the
other way the moment equations have no admissible solution. sfm() warns,
reports sigu = 0 with the whole residual variance assigned to sigv, and
returns no efficiency predictions -- to be read as no evidence of
inefficiency in these data, not as an estimate of zero.
Standard errors: the OLS slope standard errors are reported and are valid as
such. The scale parameters and the corrected intercept carry NA by
default, since neither has a closed-form standard error here and the OLS
intercept standard error would be wrong (it knows nothing about the sampling
error of a third-moment estimate). Set cols_boot for a nonparametric
bootstrap covering every parameter, with rand.cols to make it
reproducible.
sfm(model_name = "tHN") -- Student's t--half-normal. Heavy-tailed
noise (v ~ sigma_v * t_nu) with a conventional half-normal inefficiency
term (u ~ |N(0, sigma_u^2)|), drawn independently.
This is not THT. In THT (Tancredi 2002) both components come from one
shared scale mixture, so both are t with the same degrees of freedom and the
composed error is a closed-form skew-t. In tHN the tails differ, there is no
closed form, and the density is the convolution
f(e) = integral_0^Inf f_v(e+u) f_u(u) du, evaluated by Gauss-Legendre
quadrature. tHN is therefore the natural parametric comparison for the
density-power robust estimators (robust = "mlqe"/"psi"/"mdpd"), which THT
cannot be, because its inefficiency term is heavy-tailed too. Parameters are
reported as (sigv, sigu, nu) -- the conventional order, not THT's
inverted one. Returns exp_u_hat and u_hat by a Bayes rule over the same
quadrature nodes.
Two documented properties, both surfaced rather than hidden. The degrees of
freedom are weakly identified: on data simulated from the model at
n = 1000 with a true nu = 5, the profile log-likelihood moves only about
0.24 across nu from 10 to 100, and peaks near 20. Profile over a grid of
fixed nu rather than reporting one selected value. Because of that flat
ridge every tHN fit runs from several widely separated starts, keeps the
best, and reports the outcome in thn_starts, warning when the starts reach
different optima. sigma_u can collapse onto zero on real data, the heavy
noise tail absorbing the whole one-sided component and leaving mean predicted
efficiency near one; sfm() warns and sets thn_sigma_u_at_bound rather than
bounding sigma_u away from zero, because the collapse is a property of the
model and is the thing a user needs to see.
The quadrature node count scales with sigma_u/sigma_v and is not fixed. A
fixed 96-node rule is accurate near lambda = 3 but carries 4% relative
error at lambda = 20 and 60% at lambda = 62 -- and the model does reach
that region. Note that integrating the density to 1 does not detect this: the
error redistributes across the support and integrates away, so total mass
still reads 1.000 while the density is 40% wrong pointwise.
data_gen_cs() gained y_pcs_thn (with v_thn and u_thn), the
matching generator for tHN.
psfm(model_name = "TFE") -- Greene (2005) true fixed effects. The
composed-error likelihood with one intercept per individual, estimated as a
profile likelihood in (lambda, sigma, beta) with the firm effects
concentrated out. Reports the same parameter layout as "TFE_WMLE", plus
r_hat_m (the maximum-likelihood firm effects), exp_u_hat and u_hat.
Note that this likelihood always has a supremum on the sigma_v = 0
boundary, because the individual effects are unrestricted. The new argument
tfe_lambda_max (default 100) bounds the search accordingly, and a fit that
pins at the bound warns. See ?psfm for the details; this is a property of
the estimator and is one of the motivations for "TFE_WMLE".
psfm(model_name = "K1990") and "K1990modified" -- Kumbhakar (1990)
time-varying inefficiency, B_it = (1 + exp(bt + ct^2))^-1 and
B_it = 1 + d(t - T_i) + e(t - T_i)^2. These share one likelihood with
"PL80" and "BC92", differing only in B_it. K1990's b and c are
weakly identified, so the fitted B_it path is more interpretable than
either coefficient on its own.
Four new inefficiency distributions in sfm(): "NU" (normal-uniform),
"NGE" (normal-generalized exponential), "NLN" (normal-lognormal) and
"NW" (normal-Weibull). The last two are estimated by simulated maximum
likelihood over Halton draws.
predict(), fitted() and residuals() methods for class "sfareg",
alongside the existing coef(), vcov(), logLik() and nobs().
predict() accepts newdata.
psfm() now accepts an ordinary data.frame (or tibble/data.table) as well
as a plm::pdata.frame; the panel index is constructed internally from
individual and the new time argument. Previously a plain data frame
failed with an uninformative "empty model" error.
psfm(collinear_action =) controls what happens when the
between-individual design used to build starting values is rank deficient
-- the situation created by, for example, time dummies, which are estimable
in a pooled specification but collapse onto the intercept once averaged
within each unit. "start_only" (default) keeps the requested model and
drops the offending columns from the starting-value regression only;
"error" stops and names them; "warn_drop" removes them from the
estimated model. Previously this surfaced as an opaque
solve(crossprod(ZBeta)) LAPACK error inside plm.
sfm(robust =) selects a divergence-based robust estimator -- "mlqe",
"psi" or "mdpd" -- with sandwich standard errors. Currently implemented
for model_name = "NHN"; other models error rather than silently ignoring
the argument.
sfm() gained use.nlminb and use.bobyqa (both "auto" by default) and
maxit.nlminb, for control over the optimizer stack described below.
nlminb now leads the optimizer stack for the models where it helps,
ahead of the derivative-free stages, with an analytic gradient supplied for
NHN. Model-by-model defaults are chosen automatically: NHN, NE, NTN
and NU use it, while NR and NGE (where it degraded the fit) do not.
Typical cross-sectional fits are several times faster with unchanged
estimates.
An nlminb stage was added to the PL80/BC92/K1990/K1990modified
branch, worth up to +25 log-likelihood on the seven-parameter models, where
psfm()'s low default maxit.bobyqa was stopping short of convergence.
BC92 now matches or beats frontier::sfa() on parameter, variance and
efficiency accuracy.
PL80 and BC92 are now estimated by a native maximum-likelihood
implementation instead of wrapping frontier::sfa(). Verified against
frontier (coefficients, log-likelihood and predicted efficiencies) on
balanced and unbalanced panels and on both production and cost
specifications before the dependency was removed.
Efficiency prediction (exp_u_hat) is now available for NE and NTN,
which previously returned none.
The vignette still described "PL80" and "BC92" as wrappers around
frontier::sfa(). They have been native closed-form maximum-likelihood
implementations since 1.1.3, the release that dropped frontier as a
dependency -- so for four releases vignette("intro_to_psfm") advertised a
dependency on a competing package that the package had deliberately removed,
and that appears in neither Imports nor Suggests. ?psfm had it right
throughout ("natively estimated ... verified against frontier::sfa()
before that dependency was removed"); only the vignette was stale.
Two related corrections in the same file. The psfm_bootstrap() section gave
the wrong reason for excluding PL80/BC92 -- it is that they do not expose
the $U/$H efficiency structure the function reads, not that they wrap
anything -- and its list of supported models omitted TFE_WMLE and
GTRE_FML, the latter being what psfm(model_name = "GTRE") returns under
its default estimator. And the model table, which says "eleven" choices,
now says twenty-one and points at ?psfm rather than implying the table is
exhaustive.
psfm_bootstrap() refused the estimator psfm() reaches by default.
Since 1.1.4, psfm(model_name = "GTRE") defaults to estimator = "fiml"
and returns an object whose $model_name is "GTRE_FML".
psfm_bootstrap() checked that name against a list that did not contain it,
so the natural two-line sequence -- fit a GTRE model, then bootstrap it --
stopped with "psfm_bootstrap() does not support model_name = 'GTRE_FML'",
and the error message's own explanation named only the moment-based and
wrapped models, none of which applied. Reaching the bootstrap at all
required knowing to pass estimator = "sml" to the fit.
"GTRE_FML" is now supported. The data-generating process is the one the
"GTRE" arm already simulated -- same four-component model, different
estimator -- so only the parameter layout differs: FIML reports the four raw
scales (sigr, sigv, sigh, sigu) after the frontier block rather than
"GTRE"'s lambda/sigma reparameterization. $H and the boundary
warning come with it.
psfm(OPG_calc = TRUE) returned NA for every OPG and sandwich standard
error, and wrote a variable into the global environment. The OPG "meat"
matrix was stored with a superassignment, OPG_meat <<- crossprod(score_mat),
on the assumption that the surrounding tryCatch({...}) introduced a scope to
escape from. It does not -- tryCatch() evaluates its expression in the
calling frame -- so <<- began its search one frame further out, skipped the
local OPG_meat <- NULL binding entirely, and assigned into globalenv().
The immediately following solve(OPG_meat) therefore still saw NULL and
failed, as did the MASS::ginv() fallback, so the OPG standard errors were
always NA; the sandwich errors, which reuse the same matrix, were NA with
them. Both failures were reported through the existing handlers as
"OPG matrix singular, using pseudoinverse" and "OPG matrix is singular",
which pointed at a rank problem in the data rather than at the scoping bug.
Changed to a plain <-. On a fixed-seed GTRE_Z fit the parameter estimates
and Hessian standard errors are bit-identical before and after, the OPG and
sandwich errors change from NA to finite values, the two warnings stop
firing, and OPG_meat no longer appears in the user's workspace.
psfm_bootstrap() failed on Windows whenever sfa was not installed in a
default library. The function distributes work over PSOCK cluster workers,
which start as fresh R sessions holding the default library path rather than
the parent session's. Where the parent had found sfa somewhere else -- a
project library, renv/packrat, a user-set R_LIBS_USER, or the temporary
sfa.Rcheck tree that R CMD check installs into -- the workers' library()
call could not see it, and the bootstrap died in
parallel:::checkForRemoteErrors() with there is no package called 'sfa'.
On Unix the workers generally inherit R_LIBS from the parent's environment,
which hid the bug; on Windows they do not. Where the workers instead found a
different, older copy of sfa in a default library, they silently ran the
bootstrap against that version rather than the one the user had loaded. The
parent's .libPaths() is now pushed to the workers before any package is
loaded, so both sessions resolve every package identically.
stats::dlnorm was used without being imported. The rewritten "NLN"
likelihood and its efficiency block call dlnorm(), but NAMESPACE did not
import it, so the call resolved only via the search path rather than the
package namespace. R CMD check --as-cran reported it as an undefined global.
Now imported explicitly.
zsfm()'s efficiency predictor used a different mixing probability from
the likelihood it maximised. "ZISF"'s likelihood sets
prob = exp(-abs(gamma)), which makes it exactly symmetric in gamma: +g
and -g fit identically and the optimizer may return either. The JLMS block
then used exp(-gamma), so a negative estimate produced a mixing
"probability" above 1 and silently invalid post.prob and jlms. This was
reachable from an ordinary starting value, not a pathology: seeding at the
negated estimate returns gamma = -0.3015 with an identical log-likelihood,
where the old code computed prob = 1.3519. Both places now use
exp(-abs(gamma)).
zsfm()'s two-component mixture is now formed on the log scale. It built
prob*exp(f1) + (1-prob)*exp(f2) and took log(f + 1e-10). Both terms
underflow to zero when an observation is unlikely under either regime, and
the 1e-10 guard then floors that observation's contribution at -23.03 --
which does not merely protect the logarithm, it makes the objective flat
across the whole region beyond the floor, exactly where the optimizer needs a
gradient. A new internal helper, .log_add2(), computes
log(exp(a) + exp(b)) without leaving the log scale, and post.prob is now
a ratio taken in logs. Fitted coefficients are unchanged to within optimizer
path noise (worst discrepancy 1.9e-3 over eight seeds, log-likelihoods
agreeing to 1e-6).
zsfm(logit = TRUE) now uses plogis() instead of
exp(eta)/(1 + exp(eta)), which overflows to Inf/Inf = NaN once the linear
predictor passes about 710 -- a value the optimizer can reach while
searching. Identical wherever the old form was finite.
zsfm()'s efficiency block now branches on model_name rather than on
is.na(n_z_vars). The likelihood already branched on the model, and the
parameter layout is a property of the model; keying the predictor off a
different condition meant the two could disagree, reading "ZISF"'s
parameters under "ZISF_Z"'s layout if n_z_vars ever arrived as 0
rather than NA.
Note left in place, not changed: the logit = FALSE branch computes
pnorm(eta)/(1 + pnorm(eta)), which is bounded above by 0.5 and is not the
probit link. Correcting it would change results for that option, which is a
modelling decision rather than a cleanup.
sfm(model_name = "NLN") was integrating a spike. Its simulated
likelihood averaged the normal kernel over lognormal draws, but the kernel is
only sigma_v wide in u while the lognormal spreads over decades, so
nearly every draw landed where the kernel is numerically zero and a handful
carried the whole integral. Measured against a reference verified two ways
(adaptive quadrature and a 200,000-point Simpson rule, agreeing to 5e-9), the
simulated log-likelihood at the true parameters was 226.8 units low at
n = 3000 under the default draw count -- which is precisely the unexplained
228-unit gap recorded against this model in the convergence notes. It was
simulation error, not a defect in the likelihood.
sfm() now substitutes u = sigma_v*t - e, turning the integral into a
standard-normal expectation of the smooth lognormal density truncated to
u > 0. The error at the same draw count is then 0.12, and 0.86 at Nsim = 50. Raising Nsim was not an alternative: the old error per observation was
about -0.076 at both n = 1000 and n = 3000 under the 8*sqrt(n) rule, so the
total bias grew linearly in n at the same rate as the log-likelihood itself,
and closing it needed Nsim proportional to n -- quadratic work per
evaluation. "NW", which uses the same machinery but a far lighter-tailed
inefficiency term, is unchanged.
sfm(model_name = "NR") was started from a flat guess it could not
recover from, and was documented as the wrong model. Two separate
problems, one in the estimator and one in everything written about it.
The estimator: start_cs() hard-codes sigma_u = sigma_v = 0.1 for the
cross-sectional models, and from there "NR" converged to a point with a
worse log-likelihood than the true parameter vector in 9 of 14
replications at n = 4000. It failed in two modes -- a hard collapse to
sigma_u = 1e-7 with sigma_v inflated to absorb the spread (a ~50
log-likelihood deficit), and a partial stall at sigma_u around 0.85 with
the intercept ~0.4 too low. "NR" is now started by inverting the Rayleigh
moment equations instead, which reaches the same optimum a truth-seeded run
finds in all 14. Because the Rayleigh skewness is a constant, the third
central moment of the residuals identifies Var(u) outright. Wrongly skewed
residuals leave the moment equations with no admissible solution, and the
old flat start is then used unchanged.
The documentation: "NR" had been described in this package as "an
alternative closed-form derivation of the same normal/half-normal composed
error" as "NHN". That is wrong. "NR" is normal-Rayleigh -- its coded
density reproduces a numerical normal-Rayleigh convolution to 1e-8 and
misses the half-normal one by 86%. The two are different families, not
reparameterizations: the standardized skewness the inefficiency contributes
is a different constant in each (-0.631 against -0.995), and no
transformation of a two-scale family moves a standardized moment. The
likelihood itself was correct throughout and is unchanged.
data_gen_cs() gains u_r and y_pcs_r to test "NR" against its own
data-generating process; it had previously been tested against y_pcs,
which it cannot fit. The new columns are appended at the end of the
function, so every existing column is bit-for-bit unchanged. sigma_u is on
the convention E[u^2] = sigma_u^2, matching "NHN", so the Rayleigh scale
is sigma_u/sqrt(2).
sfm(model_name = "NG") and "NNAK" read their frontier coefficients
from the wrong position in the parameter vector. The likelihood closure
slices the coefficients out by a fixed offset -- x[3:(n_x+2)] for models
with two leading scale parameters, x[4:(n_x+3)] for three. NG and NNAK
carry three (sigv, sigu, mu) but sat in the two-parameter group, so
the slice took the right number of coefficients starting one slot too
early.
The consequences were severe and entirely silent. mu was used
simultaneously as the gamma (or Nakagami) shape and as the intercept
coefficient, so it could not move freely -- which is why the shape appeared
never to leave its starting value. Every remaining slope was shifted one
place, and the last coefficient never entered the likelihood at all, so it
simply kept whatever starting value it was given. The efficiency block used
the correct offset throughout, so the two halves of the model disagreed
about which number meant what.
On the package's own test DGP at n = 4000, NG returned sigma_u = 0 and
stopped 710 log-likelihood units below the true parameter vector. Across
15 fits spanning n = 1000 to 5000, none reached the truth's likelihood.
After the fix, all 15 do, no fit collapses, and mean RMSE against the truth
falls from 0.60 to 0.27 and declines with n. NNAK improves on the same
fix -- 6 of 7 fits now beat the truth and standard errors are finite, where
previously the Hessian was singular in every replication -- but it still
produces occasional failures and is not yet considered repaired.
The NG density itself was never wrong: it agrees with numerical convolution of the normal and gamma densities to 7e-14 across the whole residual range, and that check is now a test.
sfm(model_name = "NG") also started in the wrong place. start_cs()
hard-codes sigma_u = sigma_v = 0.1 and the NG start pinned the shape at 1,
so the search began from E[u] = 0.1 against a true 1, with the intercept
at the raw OLS value -- itself E[u] below the frontier. NG now builds
candidate starts from the residual moments and sweeps the shape along the
E[u] = mu*sigma_u ridge, which is the direction the data leave weakly
determined (Ritter and Simar, 1997), polishing the most promising few before
choosing. The search is reported in $ng_starts.
psfm(model_name = "GTRE_FML") started its search from the wrong
intercept. It seeded beta_0 at the raw panel-regression intercept, while
"GTRE" and "TRE" seed theirs at that intercept plus E[u] + E[h]. In a
composed-error model the regression intercept sits below the frontier by
(sigma_u + sigma_h) * sqrt(2/pi) -- 1.12 at the package's own test DGP --
so the FIML search began a full unit low, in exactly the direction of the
sigma_h = 0 boundary optimum where the model collapses to "TRE" and the
intercept absorbs the missing E[h].
Fits that fell in were not merely imprecise: on one of six test draws the
reported solution had a lower log-likelihood than the true parameter
vector (-5267.87 against -5267.80), with sigma_h = 2e-16, beta_0 = 0.155
against a true 0.5, and sigma_r inflated to 0.307 against a true 0.2 --
a local optimum, not an estimate.
"GTRE_FML" now also builds a second candidate start from the two-step
moment estimator (the one "GTRE_SEQ2" reports), evaluates the likelihood
at both and begins from the better, as Colombi (2010) and Colombi, Martini
and Vittadini (2011) recommend for this likelihood. Across the six test
draws, boundary collapses go from one to none, mean RMSE against the truth
falls from 0.049 to 0.023, and all six fits now reach a higher likelihood
than the truth. The chosen start is reported in start_search.
sfm(model_name = "THT") used the wrong likelihood. The skew-t density
of Tancredi (2002, eq. 4) has scale omega = sqrt(sigma_v^2 + sigma_u^2),
but the implementation evaluated the Student-t factor at the raw residual
and omitted the 1/omega Jacobian, which pins the scale at 1. Because
2*f(e)*G(w(e)) is a valid density for any symmetric f and odd w
(Azzalini's lemma), the wrong version still integrated to 1 and still
produced plausible fits rather than an obvious failure -- but sigma_u and
sigma_v were then identified only through the skewing term, and the
degrees of freedom a had to absorb the scale mismatch. Fitted a and
sigma_v were consequently inconsistent. Fixed; THT now reproduces
sn::dst() to machine precision.
sfm(model_name = "THT") now reports efficiency. It previously returned
no efficiency prediction at all. It now returns exp_u_hat = E[exp(-u)|e],
u_hat = E[u|e], and sd_exp_u_hat, following Tancredi (2002, section 2.2).
The conditional density in that paper's eq. (7) is a Student-t truncated to
the non-negative half-line with df = a + 1, location -e*sigma_u^2/omega^2
and scale sqrt((a + e^2/omega^2)*sigma_v^2*sigma_u^2/(omega^2*(a+1))),
which reduces to the Jondrow et al. (1982) normal predictor as a grows.
Note the behaviour this buys: for a large positive residual the
half-normal model drives predicted efficiency to 1 with near-zero
uncertainty, whereas the skew-t model reads the point as an outlier, so
efficiency turns back down and sd_exp_u_hat widens.
sfm(model_name = "THT") degrees of freedom are better started and
bounded. The starting value was 1 -- the Cauchy case, which has neither a
mean nor a variance -- and the lower bound was 1e-7. The start is now a
moment estimate from the excess kurtosis of the OLS residuals
(a ~ 4 + 6/kurtosis, clipped to [3, 30]), and the lower bound is 2.05, so
the search stays where the composed error has both moments.
data_gen_cs() gained y_pcs_st, which is the column THT should be
tested against. The existing y_pcs_t draws its two error components as
two independent rt() variates; that shares the degrees of freedom but not
the mixing variable, so the composed error is not skew-t and THT cannot
recover its own parameters from it. y_pcs_st uses the single common
Gamma(a/2, a/2) mixing of Tancredi eq. (5), with lam_st, u_st and
v_st also returned. y_pcs_t is retained unchanged, because renumbering
the random draws would alter every column generated after it.
nobs() on an sfm(), zsfm() or ttsfm() fit returned NA when called
from inside a function, because the recorded data argument was re-evaluated
in the caller's frame rather than the one the model was fitted in. This
propagated silently to BIC(), which needs the observation count, while
AIC() kept working.
logLik() on the estimators that are not maximum-likelihood
(GTRE_SEQ1, GTRE_SEQ2, SSFE) returned an unclassed NA, which made
AIC() and BIC() return numeric(0) -- a missing value that disappeared
instead of propagating. They now return NA as documented.
psfm()'s TFE and SSFE models silently depended on the first two
columns of data being the panel index, because plm() was called without
an explicit index. Any data whose first two columns were something else
produced "empty model".
Fixed a checkSymmetricPositiveDefinite() failure in the GTRE models
caused by asymmetric dimnames on an otherwise symmetric covariance matrix.
The Halton draws used by the simulated-likelihood models were reshaped in column-major order, so each observation drew from a narrow, non-equidistributed slice of the sequence.
A singular Hessian now yields NA standard errors rather than aborting the
whole fit.
All 76 stop() and warning() calls in R/ now pass call. = FALSE; 24 of
them did not, which made the error output inconsistent between older and newer
code paths.
The seven sapply() calls in psfm.R -- six extracting ridge/method
diagnostics from the GTRE posterior solver, one building the transient
efficiency vector -- are now vapply() with explicit numeric(1) and
character(1) templates, so a change in what the solver returns fails loudly
instead of silently producing a list column.
Removed the dependency on frontier, along with eight other packages.
Imports went from 25 packages to 15.
Lowered the R requirement from R (>= 4.4.0) to R (>= 4.0.0).
Added a testthat suite covering the model branches, the numerical helpers,
the S3 methods and the data generators.