Fits a cross-sectional stochastic frontier in which the noise \(v\) and the inefficiency \(u\) are dependent, with the dependence carried by a copula. Every other model in the package assumes independence.
copsfm(formula, data,
copula = c("gaussian", "fgm", "frank",
"clayton", "clayton90", "clayton180", "clayton270",
"gumbel", "gumbel90", "gumbel180", "gumbel270",
"joe", "joe90", "joe180", "joe270"), inefdec = TRUE,
n_nodes = 128, maxit.bobyqa = 10000, maxit.psoptim = 1000,
maxit.optim = 1000, start_val = FALSE, PSopt = FALSE,
optHessian = TRUE, Method = "L-BFGS-B", verbose = FALSE,
rand.psoptim = NULL)An object of class "sfareg", with copula and copula_par
alongside the usual components. jlms is \(E[u \mid \varepsilon]\)
computed by the same quadrature the likelihood used.
Two-sided formula for the frontier, y ~ x1 + x2. A
| segment is an error: variance determinants are not yet supported
alongside a copula.
A data.frame.
The dependence family. See ‘Choosing a family’ below.
| family | parameter | independence at | dependence it can express |
"gaussian" | \(\rho \in (-1,1)\) | 0 | both signs, no tail dependence |
"fgm" | \(\theta \in [-1,1]\) | 0 | both signs but weak: Spearman \(\rho = \theta/3\) |
"frank" | \(\theta \in \mathbb{R}\) | 0 | both signs, full range, no tail dependence |
"clayton" | \(\theta > 0\) | \(\theta \to 0\) | positive, lower-tail dependence |
"gumbel" | \(\theta \ge 1\) | 1 | positive, upper-tail dependence |
"joe" | \(\theta \ge 1\) | 1 | positive, heavier upper tail than Gumbel |
Clayton, Gumbel and Joe carry only positive dependence, which is a
real restriction: nothing rules out a negative association between noise and
inefficiency, and against negatively dependent data those families can only
report the independence boundary. The 90 and 270 suffixes are
the rotations that reverse the sign, and 180 is the survival copula,
which preserves it. So "clayton270" gives lower-tail dependence with a
negative association.
TRUE for a production frontier, FALSE for cost.
Gauss--Legendre nodes for the integral over \(u\).
As in sfm.
The package's taxonomy already writes the general joint density as $$f_{V,U}(v,u) = f_V(v)\,f_U(u)\,c(F_V(v), F_U(u); \rho),$$ with every existing specification setting \(c \equiv 1\). This function relaxes exactly that and nothing else: the marginals remain normal and half-normal. Composing over \(v = \varepsilon + Su\), $$f_\varepsilon(\varepsilon) = \int_0^\infty f_V(\varepsilon + Su)\,f_U(u)\, c\!\left(F_V(\varepsilon + Su), F_U(u)\right) du,$$ a one-dimensional integral evaluated by Gauss--Legendre quadrature on \(u = t/(1-t)\) rather than by simulation.
Why n_nodes defaults to 128. Measured against the closed-form
normal/half-normal density at independence, the largest absolute error in
\(\log f\) is 1.3e-2 at 32 nodes, 2.5e-5 at 64, and 8.2e-14 at 128. Sixty-four
looks adequate and is not: 1e-5 per observation is 0.02 in a log-likelihood over
2000 points, which is the scale at which competing modes are compared.
The dependence parameter needs a large sample. This is the thing to know
before using it. The frontier slopes behave normally, but \(\rho\) does not.
Most of these families do not recover their own dependence parameter, and
copsfm warns when you pick one. A correct density does not imply an
estimable parameter. Measured on 25 samples per family, generated from that
family at \(n = 400\) and refitted with it:
| family | mean estimate (truth) | fits on the independence bound |
"frank" | 5.43 (5) | 0% |
"clayton" | 2.24 (2) | 0% |
"gumbel" | 1.29 (2) | 36% |
"joe" | 1.52 (2) | 40% |
"clayton270" | 0.23 (2) | 56% |
"gumbel90" | 1.26 (2) | 60% |
On data generated from a Gumbel copula with Spearman \(\rho = 0.685\) at \(n = 2000\), every family -- including the true one -- returns the independence boundary, and their log-likelihoods differ by less than 0.04. The likelihood is close to flat in the dependence parameter. This is a property of the model rather than of the implementation: each density is checked against the second mixed partial of its own CDF, and each sampler against the family\'s theoretical Spearman correlation. The remaining rotations were not measured, and warn that they were not, rather than implying either outcome.
Prefer "frank" or "clayton", and treat an estimate from any
other family as exploratory.
Choosing a family, and a warning about doing so. The families differ in two things that matter: whether they admit negative dependence, and where they concentrate it. Gaussian and Frank spread association evenly and have no tail dependence; Clayton concentrates it in the lower tail and Gumbel and Joe in the upper; FGM can only express weak association at all.
Against that, weigh what the next paragraphs establish: the dependence parameter is estimated very imprecisely. On a Gaussian design at \(n = 600\) its sampling standard deviation is 0.385 against a truth of 0.5. Offering fifteen families does not make that better, and it makes one thing worse -- with fifteen candidates, some family will attain the highest likelihood on independent data by chance. On a sample generated with independent \(v\) and \(u\), the fitted log-likelihoods across families span only about 0.8, and the best of them reports \(\theta = 1.83\) for Gumbel. Treat a family comparison as descriptive, not as evidence of a dependence structure, and prefer a family chosen for a reason to one chosen by its likelihood.
On a Gaussian-copula design with true \(\rho = 0.6\), \(\sigma_U = 1\), \(\sigma_V = 0.4\):
| \(n\) | \(\hat\rho\) | \(\hat\sigma_U\) | \(\hat\sigma_V\) |
| 2000 | -0.247 | 0.571 | 0.275 |
| 8000 | 0.528 | 0.921 | 0.379 |
| 20000 | 0.571 | 0.971 | 0.389 |
It is consistent -- the estimates converge on the truth -- but at \(n = 2000\) the sign can come out wrong. The cause is joint identification rather than a bad optimiser: with \(\beta\) and the scales held at their true values the likelihood peaks exactly at \(\rho = 0.6\), so \(\rho\) is well identified conditionally; it is the trade-off against \(\sigma_U\), \(\sigma_V\) and the intercept that defeats it in moderate samples, and the fitted point genuinely attains a higher likelihood than the truth.
Treat \(\hat\rho\) as informative only in large samples, and compare against
an independent fit (sfm) before reading anything into it.
And below a few hundred observations, the independent fit is often the better estimate -- even when the dependence is real. Fitting the true copula gives a nearly unbiased but very noisy \(\hat\sigma_U\); ignoring the copula gives a precise but biased one. RMSE of \(\hat\sigma_U\) (truth 1.0) on Gaussian-copula data, 25 replications:
| design | \(n\) | true copula | independence |
| \(\rho = 0.5\) | 400 | 0.431 | 0.240 |
| \(\rho = 0.5\) | 1000 | 0.151 | 0.220 |
| \(\rho = -0.5\) | 400 | 0.251 | 0.196 |
| \(\rho = -0.5\) | 1000 | 0.175 | 0.203 |
The crossover sits between \(n = 400\) and \(n = 1000\). Below it the variance
cost of estimating \(\rho\) jointly with \(\sigma_U\) exceeds the bias it
removes, so plain sfm wins on RMSE; above it copsfm wins. The
independence bias does not shrink with \(n\) (about 22% of \(\sigma_U\), since
it is misspecification rather than noise), while the copula fit's RMSE falls as it
should.
Node counts, measured in paired runs. Forty datasets sent to every node count, so the differences are the quadrature and not the draw. Against the 128-node default, the mean difference in \(\hat\rho\) is +0.0002 at 256 nodes (\(t = 1.37\)) and +0.0046 at 64 (\(t = 1.47\)) -- neither distinguishable -- but -0.0783 at 32 (\(t = -2.79\), \(p = 0.008\)), with per-dataset differences as large as 0.52. So the quadrature has converged by 64, and 32 must not be used as a fast option: it is biased, not merely noisier. For scale, \(\hat\rho\)'s own sampling standard deviation is 0.385 at \(n = 600\), which dwarfs every node-count effect here.
Families. Only densities that could be verified are offered. Each is checked in the tests two ways: it integrates to 1 over the unit square, and it returns exactly 1 at the independence parameter. Frank, Clayton and Gumbel are deliberately absent rather than transcribed without a source to check against.
Smith, M. D. (2008). Stochastic frontier models with dependent error components. The Econometrics Journal, 11(1), 172--192.
sfm for the independent case.
# \donttest{
set.seed(4)
## n_nodes = 64 rather than the default 128: the quadrature is converged by 64
## (see Details), and this is the slowest entry point in the package.
n <- 600
x1 <- rnorm(n); x2 <- rnorm(n)
z1 <- rnorm(n); z2 <- 0.6 * z1 + sqrt(1 - 0.6^2) * rnorm(n)
y <- 0.5 + 0.8 * x1 - 0.4 * x2 + 0.4 * z1 - qnorm((1 + pnorm(z2)) / 2)
dat <- data.frame(y = y, x1 = x1, x2 = x2)
fit <- copsfm(y ~ x1 + x2, data = dat, n_nodes = 64)
fit$out
# }
Run the code above in your browser using DataLab