The “Get started” article introduces the copula correction using
CopRegPG(), the original estimator of Park and Gupta
(2012). This article covers the estimators that came after it, what each
one assumes, and when you would use it instead of PG.
Why PG is not the whole story
CopRegPG() assumes the endogenous regressor is
uncorrelated with the exogenous regressors. Under that assumption the
copula term, built from the endogenous regressor’s own marginal
distribution, is not stealing any identifying variation that the
exogenous regressors would otherwise supply. When the assumption fails,
and endogenous and exogenous regressors are correlated in most real
data, PG’s estimates become unreliable. The estimators covered below
exist because that assumption often does not hold, and each one relaxes
it in a different way.
A motivating example
sim_endog, the dataset shipped with the package, has an
endogenous regressor, z_endog, that is noticeably
correlated with one of the exogenous regressors,
w_instr:
data(sim_endog)
set.seed(1)
idx <- sample(nrow(sim_endog), 300)
d <- sim_endog[idx, ]
cor(d$z_endog, d$w_instr)
#> [1] 0.5630987
cor(d$z_endog, d$x_exog)
#> [1] -0.04315322z_endog and w_instr are correlated at about
0.56 in this sample, while z_endog and x_exog
are close to uncorrelated. That correlation between z_endog
and w_instr is exactly the situation PG is not built for.
validity(), covered in the next article, tests this
directly and reports whether PG’s assumption is violated in a given
fit.
2sCOPE and IMA: transform the exogenous regressors too
CopReg2sCOPE(), from Yang, Qian and Xie (2025), relaxes
PG’s assumption by transforming the exogenous regressors to normal
scores as well, not just the endogenous one. It then regresses the
endogenous regressor’s normal score on the exogenous regressors’ normal
scores in a first stage, and keeps the residual as the copula term for
the second stage. Because that residual is orthogonal to the exogenous
regressors by construction, the correlation that broke PG’s assumption
is projected out before the copula term ever enters the structural
regression.
CopRegIMA(), from Haschka (2025), runs the same two
stages and differs from 2sCOPE in exactly one place. Its first-stage
regression has no intercept, which is what the paper’s derivation
requires. The two estimators are usually close, and the gap between them
widens the further the transformed exogenous regressors sit from mean
zero, for example with an unbalanced factor.
BMW: transform the first-stage residuals instead
CopRegBMW(), from Breitung, Mayer and Wied (2024),
inverts the order of 2sCOPE’s two operations. It regresses the raw
endogenous regressor on the raw exogenous regressors first, and only
then applies the rank transform, to the residual of that regression
rather than to the regressor itself. Because the identifying assumption
now concerns those first-stage residuals rather than the endogenous
regressor, validity() tests the residuals for this
estimator instead of testing the regressor directly. The reordering also
has a second consequence. Corollary 3.2 of the paper shows that the
textbook t statistic built from classical OLS standard errors stays
valid for testing whether the copula term’s coefficient is zero, so
summary() reports that Durbin-Hausman-Wu test alongside the
usual bootstrap one.
set.seed(1)
fit_bmw <- CopRegBMW(y ~ z_endog | x_exog + w_instr, data = d, nboots = 99)
summary(fit_bmw)
#>
#> Copula endogeneity correction: BMW (Breitung, Mayer & Wied 2024)
#>
#> Call:
#> CopRegBMW(formula = y ~ z_endog | x_exog + w_instr, data = d,
#> nboots = 99)
#>
#> Residuals of the augmented regression (u = xi - C gamma):
#> Min 1Q Median 3Q Max
#> -2.65236 -0.74424 0.03662 0.69748 2.36411
#>
#> Residuals of the structural model (xi = y - mu - P alpha - W beta):
#> Min 1Q Median 3Q Max
#> -4.48245 -1.03323 0.01079 1.03488 3.78749
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 1.08412 0.10908 9.939 <2e-16 ***
#> z_endog 0.78058 1.00395 0.778 0.437
#> x_exog 1.43928 0.09951 14.464 <2e-16 ***
#> w_instr 0.88118 0.66789 1.319 0.187
#> z_endog_cop 1.16738 0.91873 1.271 0.204
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Endogeneity: rho(P*, xi*) is the correlation between the normal score
#> of an endogenous regressor and that of the structural error, xi* = xi / sigma.
#> Estimate Std. Error z value Pr(>|z|)
#> rho(z_endog*, xi*) 0.7585 0.5619 1.35 0.177
#>
#> Durbin-Hausman-Wu test of no endogeneity, rho = 0. Under that null the
#> textbook t statistic stays valid even though the copula terms are
#> generated regressors (Breitung, Mayer & Wied 2024, Corollary 3.2), so
#> these standard errors are classical rather than bootstrapped:
#> Estimate Std. Error t value Pr(>|t|)
#> z_endog_cop 1.167 1.134 1.03 0.303
#>
#> Fit, on 295 residual degrees of freedom:
#> augmented structural
#> Residual standard error 0.9950 1.5268
#> R-squared 0.8811 0.7200
#> Adjusted R-squared 0.8795 0.7162
#> sigma above is the standard error of the structural model, the one
#> entering xi* = xi / sigma.
#> Standard errors from 99 bootstrap replicates; cdf = "rank.n1", ties = "max".
#> Pr(>|z|) in both tables: Wald test using the normal approximation,
#> z = Estimate / Std. Error, with the bootstrap standard error.
#> See confint(object, type = "percentile") for bootstrap percentile intervals.
#>
#> --- Identification diagnostics ------------------------------------
#>
#> Non-normality of the endogenous regressors (small p = non-normal, good):
#> AD AD p KS p
#> z_endog 0.4423 0.2862 0.7698
#>
#> Correlation of the copula terms with the exogenous regressors
#> (this estimator does not assume it is zero; it projects the correlation
#> out in the first stage):
#> max |corr| with p (Holm) joint R2 joint p
#> z_endog 0.004219 w_instr 1 2.258e-05 0.9967
#> full matrix in summary(object)$diagnostics$exog.correlation.matrix
#>
#> Collinearity of the copula terms (omega near 0 = weakly identified):
#> corr(P, C) omega
#> z_endog_cop 0.8252 0.002654The Durbin-Hausman-Wu test block near the middle of the
output is the one this estimator adds. It is not present in the
summary() output of any of the other estimators on this
page.
JAMS: let the copula structure vary across categories
CopRegJAMS(), from Liengaard, Becker, Bennedsen, Heiler,
Taylor and Ringle (2025), lets the copula structure itself differ across
the categories formed by the discrete exogenous regressors, controlled
by the conditional argument.
conditional = TRUE, the default, uses the joint categories
of every factor or character regressor on the exogenous side.
conditional = FALSE imposes one common copula structure
over the whole sample. A character vector naming specific variables
restricts the splitting to those, useful once the joint categories
implied by TRUE would leave too few observations per cell.
sim_endog has no factor regressor, so
conditional = TRUE and conditional = FALSE
coincide here. There is nothing to split on.
set.seed(1)
fit_jams <- CopRegJAMS(y ~ z_endog | x_exog + w_instr, data = d, nboots = 99)
coef(fit_jams)["z_endog"]
#> z_endog
#> 0.9854831To see JAMS actually split the copula structure, we need a factor. The small simulation below adds a two-level region variable:
set.seed(3)
n <- 300
region <- factor(sample(c("north", "south"), n, replace = TRUE))
w <- rnorm(n)
u <- rnorm(n)
p <- rchisq(n, df = 3) + 0.4 * u
y <- 1 + 2 * p + 0.5 * w + 0.8 * (region == "south") + u
dat_region <- data.frame(y = y, p = p, w = w, region = region)
set.seed(3)
fit_jams_region <- CopRegJAMS(y ~ p | w + region, data = dat_region, nboots = 99)
summary(fit_jams_region)
#>
#> Copula endogeneity correction: JAMS (Liengaard et al. 2025)
#>
#> Call:
#> CopRegJAMS(formula = y ~ p | w + region, data = dat_region, nboots = 99)
#>
#> Residuals of the augmented regression (u = xi - C gamma):
#> Min 1Q Median 3Q Max
#> -3.220601 -0.613459 0.009187 0.651744 3.239971
#>
#> Residuals of the structural model (xi = y - mu - P alpha - W beta):
#> Min 1Q Median 3Q Max
#> -2.69694 -0.74605 0.03236 0.76292 3.28057
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 1.30839 0.21272 6.151 7.72e-10 ***
#> p 1.85385 0.08051 23.027 < 2e-16 ***
#> w 0.60114 0.06497 9.253 < 2e-16 ***
#> regionsouth 0.78742 0.12697 6.202 5.59e-10 ***
#> p_cop.regionnorth 0.47832 0.16773 2.852 0.00435 **
#> p_cop.regionsouth 0.49044 0.18641 2.631 0.00851 **
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Endogeneity: rho(P*, xi*) is the correlation between the normal score
#> of an endogenous regressor and that of the structural error, xi* = xi / sigma.
#> Estimate Std. Error z value Pr(>|z|)
#> rho(p*, xi*) 0.4591 0.1221 3.759 0.000171 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Joint significance of the copula terms (bootstrap Wald test):
#> chi-squared = 8.601 on 2 df, p = 0.01356
#> Copula structure constant across the 2 categories of the discrete
#> regressors (bootstrap Wald test): chi-squared = 0.009933 on 1 df,
#> p = 0.9206. A small p favours the grouped estimation.
#>
#> Fit, on 294 residual degrees of freedom:
#> augmented structural
#> Residual standard error 0.9664 1.0898
#> R-squared 0.9652 0.9557
#> Adjusted R-squared 0.9646 0.9549
#> sigma above is the standard error of the structural model, the one
#> entering xi* = xi / sigma.
#> Standard errors from 99 bootstrap replicates; cdf = "ecdf.adj", ties = "max".
#> Pr(>|z|) in both tables: Wald test using the normal approximation,
#> z = Estimate / Std. Error, with the bootstrap standard error.
#> See confint(object, type = "percentile") for bootstrap percentile intervals.
#>
#> --- Identification diagnostics ------------------------------------
#>
#> Non-normality of the endogenous regressors (small p = non-normal, good):
#> AD AD p KS p
#> p 9.436 6.763e-23 7.236e-05
#>
#> Correlation of the copula terms with the exogenous regressors
#> (this estimator does not assume it is zero; it projects the correlation
#> out in the first stage):
#> max |corr| with p (Holm) joint R2 joint p
#> p 0.1124 w 0.1036 0.01264 0.1513
#> full matrix in summary(object)$diagnostics$exog.correlation.matrix
#>
#> Collinearity of the copula terms (omega near 0 = weakly identified):
#> corr(P, C) omega
#> p_cop.regionnorth 0.8904 0.2743
#> p_cop.regionsouth 0.9505 0.2023The coefficient table now carries two copula terms,
p_cop.gnorth and p_cop.gsouth in the row names
(one per level of region), and summary() adds
two bootstrap Wald tests below the endogeneity table: one for the joint
significance of the copula terms, and one asking whether they differ
across the categories at all. A small p value on the second test favors
the grouped estimation that JAMS provides over a single common
structure.
2sCOPE-np: assume the least, pay the most
CopReg2sCOPEnp(), from Hu, Qian and Xie (2025), replaces
the entire first stage with a nonparametric estimate of the conditional
distribution of the endogenous regressor given the exogenous ones, using
the kernel estimator of Li and Racine (2008). It assumes nothing about
the shape of either regressor’s marginal distribution and nothing about
the functional form linking them, which makes it the estimator with the
fewest assumptions in the package. It is also the only estimator here
that tolerates a discrete endogenous regressor, such as a count or a
binary treatment, because the kernel smooths the regressor itself rather
than requiring a continuous marginal CDF to invert.
That generality costs runtime. Bandwidth selection by
cross-validation is an O(n^2) computation, and by default it is repeated
from scratch in every bootstrap replicate, so a call with a few hundred
observations and 199 bootstrap replicates can take a long time. The
chunk below is not run in this article for that reason; set
bw.refit = FALSE if you want faster, slightly less
faithful-to-the-paper standard errors, and start with a small
nboots while you are still exploring a model.
fit_np <- CopReg2sCOPEnp(y ~ z_endog | x_exog + w_instr, data = d, nboots = 99)
summary(fit_np)Unlike every other estimator on this page,
CopReg2sCOPEnp() has neither a cdf nor a
ties argument. There is no marginal CDF to choose, because
the copula term is built from a conditional distribution instead.
The cdf argument
Every estimator except 2sCOPE-np needs an estimate of a marginal CDF
to build its copula term or terms, and the literature does not agree on
a single way to get one. endogCopula implements seven,
selected with the cdf argument:
-
"kde.silverman": the integral of a Gaussian kernel density with Silverman’s bandwidth, the estimator Park and Gupta use in the original article, and the default forCopRegPG(). -
"kde.cv": the same kernel with a cross-validated bandwidth. -
"kde.plugin": the same kernel with a plug-in bandwidth. -
"ecdf.fixed": the empirical CDF with a replaced boundary, proposed by Becker, Proksch and Ringle (2022). -
"ecdf.adj": an adjusted empirical CDF, the default forCopRegJAMS(). -
"rank.n": a rescaled empirical CDF, Rank/n with a correction, the default forCopReg2sCOPE()andCopRegIMA(). -
"rank.n1": Rank/(n + 1), the default forCopRegBMW().
Each estimator already defaults to the CDF its own paper uses, so
leaving cdf alone is a reasonable starting point.
CopRegBMW() is the one place where the choice is not really
free. Its asymptotic theory is derived specifically for
"rank.n1", and supplying anything else triggers a warning
that the reported standard errors are no longer covered by that theory.
A companion argument, ties, controls how tied values are
ranked when computing any of the rank based or empirical CDFs:
"max" (the default) is the counting function used
throughout this literature, "average" uses midranks. The
two agree once there are no ties.
copreg(): one entry point for all six
copreg() reaches every cross-sectional estimator above
through a single method argument, one of "pg",
"2scope", "ima", "bmw",
"jams" or "np". It does nothing more than
rewrite the call and hand it to the estimator’s own function, so the two
routes behave identically, warnings included.
set.seed(1)
fit_generic <- copreg(y ~ z_endog | x_exog + w_instr, data = d,
method = "2scope", nboots = 99)
identical(unname(coef(fit_generic)), unname(coef(fit_2scope)))
#> [1] TRUECopRegPANEL() and CopRegBAYES(), in the
companion packages endogCopulaPanel and
endogCopulaBayes, are deliberately not reachable through
copreg(): the panel estimator needs an index identifying
the panel structure and uses a different likelihood, and the Bayesian
estimator replaces bootstrap inference with a posterior.
The coefficients side by side
Putting five of the six cross-sectional estimators on the same data and the same endogenous coefficient shows both how much they can disagree and how PG’s broken assumption shows up in practice. 2sCOPE-np is left out of the table because its cross-validated bandwidth selection is too slow to run in an article, as noted above:
ols <- lm(y ~ z_endog + x_exog, data = d)
set.seed(1); fit_pg <- CopRegPG(y ~ z_endog | x_exog + w_instr, data = d, nboots = 99, verbose = FALSE)
tab <- data.frame(
estimator = c("OLS", "PG", "2sCOPE", "IMA", "BMW", "JAMS"),
estimate = round(c(coef(ols)["z_endog"],
coef(fit_pg)["z_endog"], coef(fit_2scope)["z_endog"],
coef(fit_ima)["z_endog"], coef(fit_bmw)["z_endog"],
coef(fit_jams)["z_endog"]), 3),
std.error = round(c(NA, fit_pg$std.error["z_endog"], fit_2scope$std.error["z_endog"],
fit_ima$std.error["z_endog"], fit_bmw$std.error["z_endog"],
fit_jams$std.error["z_endog"]), 3)
)
tab
#> estimator estimate std.error
#> 1 OLS 2.090 NA
#> 2 PG 0.391 0.676
#> 3 2sCOPE 0.767 0.562
#> 4 IMA 0.767 0.562
#> 5 BMW 0.781 1.004
#> 6 JAMS 0.985 0.550PG’s estimate is far below every other estimator’s, a direct
consequence of the correlation between z_endog and
w_instr that PG is not built to handle. The copula term
ends up absorbing some of w_instr’s effect along with the
endogeneity correction. 2sCOPE and IMA, which project that correlation
out in a first stage, land close to each other, as the theory predicts.
BMW, which transforms the first-stage residuals rather than the
regressor, lands in a similar range but with a visibly larger standard
error here, because with only one exogenous regressor really driving the
first stage its residual carries less identifying signal. JAMS, whose
copula construction differs from all of the above even before any
category splitting is considered, gives a somewhat higher estimate. None
of the five is guaranteed to recover the data generating process exactly
in a single sample of 300 observations. What the table shows is that the
choice of estimator has real consequences here, precisely because PG’s
assumption is violated, and that 2sCOPE, IMA, BMW and JAMS agree with
each other far more than any of them agrees with PG.
For the identification checks behind these numbers, including a test of whether PG’s assumption is actually violated in a given fit, see the next article, “Checking the identifying assumptions”.